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20.  The  identification  procedure  consists  of  the  sequential  application  of  a model 
determination  algorithm  and  a maximum  likelihood  identification  method . Input 
design  and  measurement  system  specification  techniques  are  combined  with  this 
application  to  form  the  integrated  parameter  identification  process. 

The  model  determination  algorithm  is  based  on  least  squares  subset  regres- 
sion techniques  .V This  algorithm  determines  which  coefficients  of  polynomial  ex- 
pansions of  the  nonlinear  aerodynamic  forces  and  moments  can  characterize  the 
stall/post-stall  flight  regime.  The  procedure  not  only  estimates  which  nonlinear 
parameters  are  significant,  but  also  gives  preliminary  numerical  estimate  error 
values  of  these  parameters.  To  further  define  the  model  structure  and  parameter 
estimates,  the  regression  results  are  used  to  initialize  a maximum  likelihood  algo- 
rithm. This  implementation  of  the  maximum  likelihood  parameter  identification 
method  has  been  previously  demonstrated  to  successfully  identify  parameters  from 
flight  data  which  is  corrupted  by  gusts  and  measurement  noise.  This  algorithm 
not  only  computes  the  parameter  and  noise  statistics  estimates,  but  also  the 
confidence  intervals  for  the  estimates  as  a function  of  angle-of-attack  and  sideslip 
angle . 

- -^This  structure  and  parameter  estimation  procedure  is  applied  to  the  response 
of  an  advanced  fighter.  Results  of  the  application  of  the  identification  procedure 
to  these  data  have  shown  that  the  model  structure  and  parameters  of  the  nonlinear 
regime  can  be  successfully  identified,  even  in  the  presence  of  high  levels  of 
process  and  measurement  noise.  It  is  further  shown  that  the  estimates  obtained 
from  a data  length  with  a given  input  can  be  used  to  successfully  predict  the 
response  of  the  aircraft  to  another  input  in  the  same  flight  regime. 

The  procedures  established  and  verified  are  highly  flexible  and  may  be 
readily  applied  to  modeling  other  classes  of  nonlinear  vehicles  from  test  data. 
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I . INTRODUCTION  AND  SUMMARY 


1.1  INTRODUCTION 

The  stall/post-stall/spin  high  angle-of- attack  operating  regime  of  advanced 
high  performance  aircraft  is  currently  recognized  as  a critical  area  of  the  flight 
envelope  which  requires  more  precise  definition  [1]  . The  complexity  of  the 
dynamic  and  aerodynamic  interactions  which  characterize  such  regimes  necessitates 
advanced  testing  and  analysis  methods.  One  of  these  methods  is  that  of  system 
identification  in  which  flight  data  from  the  critical  regimes  is  processed  to  isolate 
and  quantify  the  significant  aerodynamic  contributions. 

System  identification  is  emerging  as  a powerful  technology  for  exploiting 
test  data  to  significantly  improve  and  increase  fundamental  knowledge  of  aircraft 
dynamics.  Such  knowledge  is  essential  for  more  quantitative  approaches  to  handl- 
ing quality  evaluation,  aerodynamic  modeling  for  pilot  simulators,  and  verification 
of  analytical  and  wind  tunnel  estimates  of  aerodynamic  parameters.  In  addition, 
the  continuing  analytical,  computational,  and  experimental  development  of  system 
identification  promises  the  realization  of  on-line  application  which  will  reduce 
flight  test  time  and  provide  rapid  evaluation  for  aircraft  design  or  control  modi- 
fications . 

As  the  potential  of  this  method  has  been  realized,  however,  so  also  have 
numerous  problem  areas  associated  with  limitations  of  the  system  identification 
process  itself.  Basic  approaches  to  aircraft  parameter  identification,  such  as 
least  squares,  have,  under  ideal  conditions,  demonstrated  the  fundamental  capa- 
bility to  determine  estimates  of  parameters . The  deficiencies  of  such  applications 
become  apparent,  however,  when  very  high  order  systems  with  process  and 
measurement  noise  errors  are  considered.  When  estimates  are  obtained  from 
such  applications  (and  computation  divergence  may  even  prevent  such  estimates), 
they  are  biased.  To  resolve  the  deficiencies  of  such  basic  techniques,  advanced 
methods  have  evolved.  Of  these,  one  in  particular-~maximum  likelihood--has 
been  widely  accepted  as  a theoretically  and  practically  powerful  approach. 


There  remain  implementation  problems  in  the  application  of  this  versatile 
method.  The  objective  of  this  project  has  been  to  address  one  of  the  most  impor- 
tant of  these--identification  of  the  high-order  polynomial  nonlinear  aerodynamic 
coefficients  which  are  characteristic  of  the  high  angle-of-attack  stall/post-stall/ 
spin  regime.  The  central  problem  which  has  been  addressed  is  the  resolution  of 
over-parameterization  in  specifying  the  order  of  the  polynomial  to  be  identified. 

The  results  of  such  over-parameterization  are  excessive  computer  time,  possible 
computational  divergences,  and,  most  importantly,  incorrect  parameter  estimates. 

The  developments  reported  here  include  a method  for  quickly  estimating 
the  required  polynomial  model  to  avoid  over-parameterization,  and  an  advanced 
maximum  likelihood  implementation  for  identifying  the  parameters  of  that  poly- 
nomi  al . 

The  procedure  formulated  for  this  project  is  a further  extension  of  an 
integrated  parameter  identification  process  involving  model  determination,  input 
design,  and  maximum  likelihood  identification  algorithms.  Previously  developed 
elements  of  this  procedure  have  been  successfully  applied  to  flight  test  data.  The 
validation  of  this  extended  procedure  is  obtained  on  a comprehensive  digital 
simulation  [1] . The  complexity  of  the  high  angle-of-attack  flight  regime  requires 
such  "controlled  data"  for  achieving  a useful  identification  program.  This  is 
because  a fixed  set  of  flight  data  does  not  allow  the  flexibility  required  to  determine 
specific  cause  and  effect  relationships  which  may  degrade  identification  accuracy. 

1 2 PRINCIPAL  RESULT 

The  principal  result  of  this  work  is  a method  which  quantifies  the  causes 
of  high  angle-of-attack  aircraft  responses.  This  method  has  been  formulated  by 
developing  advanced  tools  for  determining  significant  aerodynamic  contributions 
to  such  responses  and  estimating  numerical  estimates  of  these  contributions. 

The  method  is  very  flexible,  and  is  designed  to  be  applicable  to  a wide  range  of 
flight  test  requirements  for  aerodynamic  coefficient  identification. 
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This  method  consists  of  the  following  algorithms: 


1 • A model  determination  procedure  using  multiple  regression 

techniques:  This  procedure  estimates  which  aerodynamic  para- 
meters adequately  describe  the  input-output  data  for  a given  response. 
This  capability  is  important  because  it  indicates  the  parameters 
which  should  be  identified,  and  thus  avoids  problems  of  attempting 
to  identify  more  parameters  than  possible  from  the  given  data.  In 
particular,  this  algorithm  will  indicate  whether  a linear  or  nonlinear 
aerodynamic  model  is  required  to  be  identified.  In  addition,  the 
procedure  provides  initial  "start-up"  estimates  for  the  maximum  like- 
lihood algorithm.  Such  preliminary  estimates  significantly  improve 
the  accuracy  of  the  maximum  likelihood  algorithm  and  reduce  compu- 
tational time  required  to  estimate  parameters. 

2 . A flexible  maximum  likelihood  algorithm  implementation:  This 

procedure  estimates  values  of  the  aerodynamic  parameters  which  are 
generally  more  accurate  than  obtainable  from  the  simpler  model 
determination  procedure . This  algorithm  further  examines  the  sig- 
nificance of  the  aerodynamic  contributions  by  an  advanced  identifi- 
ability  criterion  (e.g.  , the  rank  deficient  solution)  . This  model 
refinement  capability,  coupled  with  the  robustness  of  the  maximum 
likelihood  algorithm  to  turbulence  and  instrumentation  errors , 
results  in  highly  accurate  parameter  estimates.  This  implementation 
is  further  highly  user-oriented  with  the  capabilities  to: 

a.  Process  lateral  modes  only,  longitudinal  modes  only,  or  coupled 
longitudinal  and  lateral  motions.  Any  subset  of  equations  of 
these  modes  may  be  deleted.  These  options  allow  the  user  to 
"tailor"  the  program  to  specific  maneuvers  which  require  only 
a few  equations  to  model  the  observed  respones . This  reduces 
computation  time,  produces  more  accurate  estimates,  and 
allows  the  user  to  isolate  particular  parameter  effects . 


- p OTPMi PWPWPPWPW^^WPWPBMBIPIWWII  • -,y-  — — ..  ■ 


b.  Include  any  combinations  of  measurements  (i.e.,  measurements 

may  be  pitch,  roll,  yaw  accelerations  or  rates,  or  both).  Various 
flight  test  programs  have  different  levels  of  available  instrumen- 
tation. Alternately,  some  instruments  may  be  inoperative  over 
a given  flight  test  record.  This  program  allows  the  user  to 
specify  which  measurements  are  available.  (The  user  must, 
however,  check  that  the  instruments  are  consistent  with  the 
parameters  to  be  estimated.) 


c.  Identify  any  coefficient  of  an  arbitrary  polynomial  up  to  a maxi- 
mal order  set  by  the  user.  This  option  is  used  with  the  model 
determination  procedure  which  indicates  the  required  linear 
and  nonlinear  aerodynamic  contributions  which  should  be  iden- 
tified. For  reasons  of  parameter  estimate  accur  icy,  the  lowest 
possible  polynomial  expansion  (e.g.  , Taylor  series  expansion 
of  the  aerodynamic  forces  and  moments)  is  desirable. 


d.  Obtain  time  histories  not  only  of  the  measurements,  but  also 

of  the  parameters  and  the  estimated  covariances  of  the  parameter 
estimates  of  the  total  aerodynamic  forces  and  moments.  This 
aspect  of  the  program  indicates  the  degree  of  success  in  esti- 
mating the  parameters  of  interest.  It  is  noted  that  measurement 
time  history  "matches"  do  not  guarantee  the  "correctness"  of 
the  parameter  estimates,  and  other  tests  (described  in  this 
report)  are  necessary. 

e.  Determine  confidence  levels  for  the  individual  coefficients  of 
polynomial  representation  of  the  forces  and  moment  coefficients. 
This  option  allows  the  user  to  determine  confidence  levels  of 
various  polynomial  coefficients  (e.g.,  derivatives)  to  the  total 
aerodynamic  forces  and  moments. 
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1 . 3 ANALYTICAL  AND  COMPUTATIONAL  RESULTS 


The  development  of  this  method  has  produced  several  significant  analytical 
and  computational  results  . These  include  the  following  analytical  extensions  to 
the  application  of  the  maximum  likelihood  algorithm: 

1.  A method  for  adaptively  modifying  the  model  structure  during  the 
identification  computations  to  reduce  the  effect  of  modeling  errors 
on  parameter  estimate  accuracy. 

2.  An  extension  of  recently  developed  techniques  in  the  optimal  input 
design  for  linear  system  parameter  identification.  This  extension 
provides  a method  to  select  those  flight  test  points  which  must  be 
used  to  identify  a coefficient  which  is  nonlinear  with  angle-of- attack . 

Computational  results  have  been  directed  toward  verifying  the  identification 
method  on  the  simulation  data.  These  results  are  summarized  in  the  following 
conclusions: 


1.  A detailed  nonlinear  simulation  can  be  used  to  discover  and  resolve 
identification  problems  in  the  high  angle -of- attack  regime.  This 
simulation  is  sufficiently  detailed  so  that  no  single  model  whose 
parameters  are  to  be  identified  can  correctly  describe  the  simulated 
response  over  any  regime.  This  is  generally  a characteristic  of 
actual  aircraft  data,  particularly  in  the  nonlinear  regimes,  which 
must  be  retained  so  that  modeling  errors  are  realistically  simulated. 
The  simulation  used  for  this  project  is  documented  in  Ref.  [l]  and 
reviewed  in  Chapter  II  of  this  report. 

2.  The  model  whose  parameters  are  to  be  identified  should  be  determined 
from  the  aircraft  response  data.  The  model  determination  program  is 
a well-defined  and  repeatable  technique  for  isolating  and  estimating 
values  of  significant  linear  and  nonlinear  aircraft  model  parameters 
for  this  objective . This  procedure  is  sufficiently  robust  to  yield 
reliable  results  in  the  presence  of  high  levels  of  both  measurement 
and  process  noise. 
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The  use  of  a model  determination  "pre-identification"  filter  reduces 
problems  of  over-parameterization  (more  parameters  than  required 
to  describe  the  response)  . Alleviation  of  such  problems  produces 
more  statistically  reliable  estimates  and  eliminates  much  extra  compu- 
tation which  produces  only  marginal  improvements. 

The  combined  model  determination  and  maximum  likelihood  identifica- 
tion algorithm  is  generally  successful  in  providing  accurate  parameter 
estimates  of  significant  nonlinear  coefficients.  There  are  three 
principal  methods  to  validate  the  accuracies  of  these  model  and  para- 
meter estimates: 

a.  Measurement  time  history  match  between  actual  and  estimated 
response  from  the  identified  model.  This  is  a necessary,  but 
not  sufficient,  criterion  of  validation. 

b.  Correlation  between  estimated  coefficients  and  those  which 
are  known.  For  validation  on  simulation  data,  this  is  a para- 
mount test.  For  similar  validation  on  flight  data,  a compari- 
son with  wind  tunnel  data,  or  previous  flight  test  results, 
may  be  made,  at  least  where  reliable  previous  data  does  exist. 


c . Prediction  of  response  capabilities  (perhaps  the  most  signifi- 
cant test)  . This  criterion  is  specified  as  the  ability  of  the 
estimated  model  from  one  flight  test  (with  a given  input)  to 
predict  the  response  of  another  flight  test  at  roughly  the  same 
condition  (with  a different  input)  . 

The  integrated  model  structure  and  parameter  identification  methods 
developed  here  are  validated  by  all  three  of  these  criteria. 

5.  Where  more  accurate  estimates  of  the  parameters  are  required  than 
are  obtained  from  a calculation  of  parameter  estimates  from  data, 
input  design,  estimate  information  aggregation  and  data  length  can 
be  modified  to  provide  improved  results  on  subsequent  calculations. 
This  improvement  is  demonstrated  in  this  report. 


Input  design  strategies  for  polynomial  nonlinear  systems  may  be 
performed  as  a series  of  sequential  linear  designs  for  those  regimes 
where  linear  terms  predominate.  Such  regimes  do  exist  even  for 
high  angle-of- attack  maneuvers  and  can  be  verified  by  the  model 
determination  algorithm. 


1.4  SUMMARY  OF  REPORT 


The  techniques  discussed  in  this  report  presume  familiarity  with  advanced 
statistical  theory,  parameter  identification  methods,  and  aerodynamics.  To  dis- 
cuss the  preliminary  assumptions  used  in  this  work , Chapter  II  is  devoted  to  a 
discussion  of  the  integrated  parameter  identification  procedure,  including  a re- 
view of  the  details  of  the  simylation . This  simulation  itself  is  extensively  document- 
ed in  Ref.  [1] . 


Chapter  III  discusses  the  basic  ideas  of  the  multiple  regression  approach 
and  outlines  the  application  of  this  method  to  the  model  determination  problem. 
(Appendix  B is  adjunct  to  this  chapter  and  discusses  the  detailed  theory  of  sub- 
set regression .) 

Chapter  IV  summarizes  the  maximum  likelihood  algorithm  implemented  for 
this  development.  (Appendix  A is  complementary  to  this  chapter  and  presents 
a unified  discussion  of  maximum  likelihood  assumptions,  equations,  and  applica- 
tions . ) 


Chapter  V discusses  input  design  for  linear  and  nonlinear  systems. 
(Appendix  C reviews  previous  advanced  developments  for  input  design  and  new 
application  to  nonlinear  systems.) 


Chapter  VI  is  a review  of  the  significant  results  produced  by  applying  the 
integrated  identification  approach  to  the  nonlinear  aerodynamic  phenomena  of  the 
high  angle-of-attack  stall/post- stall  regime. 

Chapter  VII  presents  conclusions  and  suggestions  for  further  development 
and  flight  test  data  acquisition . 


^ — 


II . THE  INTEGRATED  PARAMETER  IDENTIFICATION  PROCESS 
2.1  BACKGROUND 

Aircraft  parameter  identification  is.  the  process  of  extracting  numerical 
values  for  the  aerodynamic  stability  and  control  coefficients  from  a set  of  flight 
test  data  (e.g. , a time  history  of  the  flight  control  inputs  and  the  resulting  aircraft 
response  variables).  Although  the  fundamental  theoretical  basis  of  identification 
has  existed  for  over  75  years,  practical  application  of  this  technology  to  aircraft 
flight  testing  has  been  attempted  only  over  the  last  three  decades.  Most  of  this 
application  has  been  limited  only  to  identification  of  low-order  linear  aircraft 
mathematical  models,  at  low  angles-ofrattack  and  Mach  number. 

There  are  three  principal  elements  of  aircraft  parameter  identification: 

(1)  the  data  processing  algorithms  (identification  method),  (2)  the  aircraft  instru- 
mentation, and  (3)  the  flight  control  inputs  (Figure  2.1).  These  elements  are 
highly  interactive. 

In  order  to  develop  a comprehensive  identification  technology  for  application 
to  the  high  angle-of- attack  flight  regime,  a method  of  approach  has  been  imple- 
mented which  integrates  these  interactive  elements.  As  diagrammed  in  Figure 
2.2,  this  method  of  approach  is  composed  of  two  basic  phases.  Phase  1 is  a detailed 
simulation  of  a high  performance  aircraft  (as  described  in  ref.  [1]).  Phase  2 
is  the  application  of  advanced  parameter  identification  methods  to  determine  the 
effects  of  inputs,  measurement  error,  and  the  algorithms  in  processing  high 
angle-of-attack  data. 
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2.2  REVIEW  OF  HIGH  ANGLE-OF-ATTACK  CHARACTERISTICS  AND  HIGH 
ANGLE-OF-ATTACK  DATA  SIMULATION 

An  extensive  discussion  of  the  dynamic  and  aerodynamic  characteristics 
of  the  high  angle-of-attack  flight  regime  for  high  performance  aircraft  is  given 
in  ref.  [1] . The  principal  motivation  for  the  interest  in  this  regime  is  the  non- 
linear forces  and  moments  which  occur  due  to  the  stall  of  airframe  components . 
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Figure  2.1  The  Integrated  Aircraft  Identification  Process 


Typical  variations  of  such  forces  and  moments  are  shown  in  Figures  2.3 

through  2.5  . Such  static  and  dynamic  coefficient  variation  with  angle-of- attack 

play  significant  roles  in  determining  the  stability  and  control  of  aircraft  in  the 

stall/post-stall  regime.  The  nonlinear  C versus  a of  Figure  2.3  leads  to  the 

m 

well-known  "pitch-up"  phenomenon  of  early  century  series  aircraft.  It  is  still 
encountered  on  aircraft,  although  recent  designs  have  reduced  its  occurrence. 
The  nonlinear  lateral/directional  characteristics  of  Figures  2.4  and  2.5  are  indi- 
cative of  the  source  of  destabilizing  forces  which  may  cause  inadvertent  spin 
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Figure  2.2  Implementation  of  the  Design  of  an  Integrated  Parameter 

Identification  Process 
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of  high  performance  aircraft.  Figure  2.6  dramatizes  the  effect  of  such  nonlineari- 
ties on  an  F-4  aircraft. 


-10  0 10  20  30  40  50  60  70  80  90  100  110 

a He9 

Figure  2.3  Variation  of  C with  a as  Derived  from  Wind 

m 

Tunnel  Tests 


As  demonstrated  in  ref.  [l] , responses  such  as  those  of  Figure  2.6  may 
be  simulated,  at  least  qualitatively,  by  digital  simulation.  In  particular, 
responses  such  as  pitch  up  and  yaw  departure  (due  to  unstable  C ) are  strongly 

nP 

dependent  on  wind  tunnel  measurable  static  nonlinearities  and  these  critical 
regions  can  be  predicted  for  flight  conditions  simulated  by  the  wind  tunnel  configur- 
ation. Unfortunately,  many  flight  conditions  may  only  be  approximated  in  wind 
tunnels  and,  in  fact,  some  aerodynamic  (and  dynamic)  nonlinearities  may  not 
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Figure  2.6  Rolling  Departure  and  Recovery  Rolls 
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be  detected  from  such  tests  due  to  scaling  or  tunnel  model  support  limitations. 

In  addition,  high  performance  maneuvers  will  induce  complicated  interactions 
between  the  pilot,  the  flow  fields,  and  the  aircraft  which  may  accentuate  destabiliz- 
ing forces . 

The  simulation  described  in  ref.  fl]  was  designed  as  a generator  of 
data  which  is  characteristic  of  that  from  the  responses  of  high  performance 
aircraft.  The  characteristics  of  this  simulation  are  summarized  in  Table  2.1. 


TABLE  2 . 1 

Characteristics  of  Aircraft  Data  Simulation 


Dynamic 

Equations 

Six  degree-of-freedom  coupled  nonlinear  equations  including 
engine  gyroscopic  effects,  non-symmetric  mass  loading  effects, 
and  asymmetric  thrust.  Euler  angles  obtained  from  integrating 
direction  cosine  matrix,  [ref.  3,4] 

Aerodynamic 

Data 

Tabular  data  from  wind  tunnel  tests  of  a model  of  subject  air- 
craft [ref  5-6] . Data  consisted  of  complete  static  and  dynamic 
force  and  moment  coefficients  for  an gle-of- attack  range 
(-10  < a < 110°) , and  sideslip  (-40  < p < 40)  . Lagrangian 
formulae  used  for  interpolation  when  integrating  equations. 

Measurement 

System 

Detailed  model  of  position,  rate,  and  acceleration  gyros,  angle- 
of-attack  and  sideslip  vanes,  accelerometers,  and  pitot  tube. 

All  measurements  subjected  to  random  errors--bias  error, 
scale  factor  error,  white  noise,  [ref.  7] 

Control 

System 

Detailed  model  of  the  SAS  of  subject  aircraft,  [ref.  8]  Autopilot 
used  to  emulate  pilot  responses. 

Random 
Disturb  ances 

Random  angle-of-attack  and  sideslip  disturbances  to  simulate 
gusts.  Random  variation  of  normal  lift  coefficient  to  simulate 
buffet. 
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The  simulation  was  subsequently  compared  with  reported  results  of  the 
stall/spin  tests  of  anF-4  aircraft  [2].  Responses  such  as  pitch-up,  yaw 
departure,  and  spin  were  obtained  with  the  simulation.  In  addition,  it  was  found 
that  wing  rock  would  be  simulated,  and  that  the  roll,  sideslip,  and  aileron  deflec- 
tion amplitudes  and  angle -of- attack  occurrence  were  comparable  with  that  recorded 
in  the  flight  test.  However,  the  frequency  of  the  simulation  was  7.5  seconds 
compared  to  4.5  seconds  for  the  actual  test.  This  frequency  discrepancy  was 
ascribed  to  uncertainties  in  the  parameters  of  the  aircraft  being  tested,  the  pilot 
responses  of  the  tests,  and  the  differences  between  the  wind  tunnel  and  actual 
values  of  dynamic  derivatives  such  as  Cg  . Based  on  the  correlation  of  the  simulated 

P 

responses  with  this  and  other  actual  responses,  the  simulation  was  deemed  accept- 
able as  a generator  of  typical  data.  Further  extensive  correlation  of  the  simulation 
was  beyond  the  scope  of  the  program  development. 

2.3  SELECTION  OF  MODELS  FOR  PARAMETER  IDENTIFICATION 

2.3.1  Simulation  Versus  Identification  Models 


As  discussed  in  ref.  [1]  and  summarized  in  Section  2. 1 of  this  report, 
the  simulation  used  to  generate  data  is  an  extensive  and  detailed  representation 
of  the  aircraft  and  data  error  sources.  Not  only  the  dynamic  equations,  but 
also  the  measurements  are  nonlinear . 

In  particular , it  is  noted  that  the  aerodynamic  tables  of  the  simulation  are 
non-analytic  functions  of  angle-of-attack , a,  and  sideslip  angle,  p.  These  data 
are  not  globally  fitted  with  polynomials  or  other  analytic  functions  of  a or  p in 
the  simulation.  There  are  two  reasons  for  retaining  the  "table  look-up"  character 
of  these  data.  First,  such  curve  fitting  could  introduce  additional  error  into  the 
simulation  as  well  as  increase  the  development  time  for  the  simulation.  Secondly, 
it  was  anticipated  that  analytic  polynomials  would  be  the  basic  identification  model, 
and  it  was  desired  not  to  identify  the  parameters  of  the  model  using  the  same 
function  which  generated  the  data. 
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This  modeling  aspect  of  the  efforts  reported  here  is  unique  in  the  application 
of  identification  algorithms  to  simulated  data.  Typical  approaches  assume  a polyno- 
mial form  for  nonlinear  parameter  variation,  generate  the  time  history  data  with 
this  model,  and  then  attempt  to  reconstruct  the  polynomial  from  the  data  [9-11]  . 
Such  an  approach  is  certainly  useful  for  testing  of  a program,  but  may  tend  to 
place  more  confidence  in  the  effectiveness  of  the  algorithm  than  is  justified. 

This  follows  because  identification  of  the  same  functional  form  which  generates 
the  data  disregards  the  modeling  error  which  may  occur  by  approximating  the 
actual  function  (unknown)  by  an  assumed  or  a priori  function.  The  approach 
used  in  this  work  is  to  generate  the  data  by  an  aerodynamic  model  which  is  of 
higher  order  than  the  model  which  is  identified.  Hence,  the  integrated  para- 
meter identification  process  of  Figure  2.2  contains  the  element  of  modeling  error 
effects  explicitly. 

2.3.2  Considerations  in  Selecting  Identification  Models 

2 . 3 . 2 . 1 The  Polynomial  Assumption 

The  use  of  polynomials  as  the  basis  of  the  identification  model  is  itself 
an  assumption  about  the  physics  of  the  aircraft  aerodynamics.  The  assumption 
is  historically  based  on  the  dependence  between  force  and  moment  coefficients 
and  independent  variables  (a  or  p,  for  example)  which  is  observed  in  wind  tunnels 
(see  Figures  2.3  - 2.5) . Recent  work  in  England  [12]  has  demonstrated  the 
validity  of  such  approximations  with  actual  aircraft  responses  in  the  subsonic 
regime.  Mathematically,  polynomial  representations  result  from  series  expansions 
about  some  reference  point  (e.g.  , trim)  and  herein  lies  the  inherent  assumption 
of  the  polynomial  approximation.  Specifically,  it  is  assumed  that  there  are  contin- 
uous derivatives  or  rates  of  changes  of  derivatives  (to  an  arbitrarily  high  order)  . 

In  the  transonic  regime  especially,  such  continuity  may  be  violated  forcing  dis- 
continuous representations.  For  this  work,  continuity  is  assumed. 

2. 3. 2. 2 Axis  System  Selection  Considerations 

The  simulation  used  to  generate  data  is  written  in  the  body  axis  coordinate 
system.  The  available  wind  tunnel  data  was  consistent  with  this  representation. 
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However,  it  must  be  noted  that  this  is  not  the  only  axis  system  which  can  be 
used,  nor  is  it  used  by  all  flight  test  agencies  [13] . 


Alternate  axes  systems  are  the  flight  stability  axis  system,  the  wind  tunnel 
axis  system,  and  the  wind  axis  system  [3,4].  These  axes  systems  are  defined 
relative  to  the  aircraft  velocity  vector.  The  body  axis  system,  the  principal 
axis  system,  and  the  so-called  instrumentation  axis  system  ore  fixed  in  the  aircraft. 
The  stability  axis  is  a hybrid  velocity  vector- aircraft  fixed  axis  system  in  that 
the  steady  state  relative  velocity  vector  defines  a body  fixed  system  (x  axis  along 
V)  for  perturbation  analysis  following  a disturbance.  A body  fixed  axis  system 
is  usually  used  for  the  moment  equations  (i.e.  , p,  q,  r)  because  the  inertias 
are  then  constant  as  would  not  be  the  case  for  the  velocity  vector  systems.  The 
translational  equations , however , may  be  equivalently  represented  in  either 
a velocity  vector  system  or  a body  fixed  system.  In  particular,  the  body  axis 
translational  states  (u,  v,  w)  are  related  to  the  wind  axis  translational  states 
(a,  p,  V)  by  the  definitions 


a = tan 


-i  _w 
u 


ft  . -1  V 

P = sin  y- 


(2.1) 

(2.2) 


For  small  angle-of- attack , these  two  systems  are  simply  related  by  neglect- 
ing higher-order  terms  so  that 


I 


a ~ 

P ‘ 


w/u 


v/V 


(2.3) 

(2.4) 


and  the  two  systems  are  equivalent,  the  states  of  the  two  systems  being  (a,  p, 
V)  or  (u,  v,  w)  . The  (a,  p,  V)  system  is  advantageous  because  all  the  states 
are  measured  directly  (e.g.  , angle-of-attack  vane,  sideslip  vane,  pitot  tube). 
The  advantage  is  that  the  identification  model  is  simply  formulated  since  the 
states  are  directly  measured.  On  the  other  hand,  measurements  of  u , v and 
w are  difficult,  at  best,  and  must  be  found  by  integrating  accelerometer  inputs 


or  by  using  Eqs . (2.1)  - (2.2)  to  find  u,  v,  and  w (since  V2  = u2  + v2  + w2) 
given  a and  p.  These  fall  into  the  category  of  derived  measurements  of  u,  v,  and 
w from  possibly  noisy  measurements  of  a and  p, 

The  advantage  of  using  the  (a,  p,  V)  system  is  offset  by  the  complexity 
of  these  equations  for  high  angle-of-attack  (or  large  sideslip  angle).  Table 
E.  1 of  Appendix  E shows  the  (a,  p,  V)  equations  for  this  situation  (compare 
with  equations  for  u,  v,  w in  Appendix  D)  . It  is  clear  that  a signficant  increase 
in  equation  complexity  occurs  with  this  representation.  Some  of  this  complexity 
occurs  because  of  the  retention  of  body  fixed  forces  such  as  gravity  and  thrust 
in  the  body  axis  frame.  Additional  complexity  occurs  because  of  coupling  of 
the  primary  longitudinal  and  lateral  static  and  dynamic  aerodynamic  forces  (eqs. 
for  V and  p)  . An  exact  formulation  in  the  u,  v,  w system  also  causes  such  coupling, 
but  only  through  secondary  aerodynamic  terms  such  as  C , etc.  It  is  also  easily 

shown  that  the  model  of  the  accelerometer  is  increased  in  complexity  (Table  E.2). 

In  summary,  instrumentation  is  referenced  to  a body  fixed  coordinate  system, 
and  aerodynamic  forces  are  increasingly  being  referred  to  a body  fixed  frame. 

It  is  concluded  that  the  u,  v,  w equations  should  serve  as  the  basis  of  the  identifi- 
cation model  for  high  angle-of-attack  analysis.  Such  a representation  avoids 
increased  complexity  of  the  equations  due  to  transformations  from  the  body  to 
wind  axis  frames  which  is  explicit  in  the  a,  p,  V system. 

It  must  be  noted,  however,  that  there  is  no  essential  difference  between  the 
linear  forms  of  the  two  representations  at  low  angles-of-attack . The  desire  to  avoid 
complexities  of  explicit  coordinate  transformations  for  high  angle-of-attack  is  based 
on  the  need  for  an  efficient  computational  structure  which  avoids  excessive  multi- 
plications . 

2.4  IDENTIFICATION  OF  MODELS 

Having  determined  an  analytical  a priori  form  for  the  aerodynamic  forces 
and  moments,  and  selected  the  principal  axis  systems,  the  framework  is  established 
for  estimating  the  structure  and  parameters  appropriate  to  a given  data  length. 

This  system  estimation  is  the  central  objective  of  this  work.  It  is  within  the 
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a priori  structure  that  a systematic  procedure  is  applied  to  selection  of  the  identifi- 
cation model  (structure  and  parameters),  parameter  estimation,  and  model  verifica- 


For  this  work,  careful  distinction  between  "models"  is  required.  To  sum- 
marize these,  the  following  clarification  is  made. 

• Simulation  Model:  This  is  the  complete  nonlinear  dynamic  and  aerodynamic 
simulation  discussed  in  Section  2.1.  This  model  generates  the  data  con- 
sidered representative  of  actual  flight  test  data.  Changes  to  this  model 

are  made  to  provide  experimental  degrees-of-freedom  for  investigating  the 
performance  of  the  subsequent  data  processing  under  various  conditions. 

• Polynomial  Model:  This  is  basically  the  same  system  dynamic  model  as  is 
used  in  the  simulation,  with  high  order  polynomials  representing  the  aero- 
dynamic force  and  moment  expansions  (as  opposed  to  the  discrete  wind 
tunnel  data  of  the  simulation  model)  . This  will  be  discussed  in  detail  in 
Chapter  III. 

• Identification  Model:  A lover-order  poiynomial  model,  this  representation 

is  the  model  which  is  actually  identified.  The  parameters  so  included  are 
a subset  of  the  parameters  of  the  more  complete  polynomial  model.  These 
subset  parameters  are  those  which  are  most  significant  in  determining 

the  response  of  the  system.  This  identification  model  parameter  solution 
technique  will  also  be  detailed  in  Chapter  III. 

The  process  of  arriving  at  the  "identification  model"  from  the  data  (generated 
by  the  simulation  or  by  flight  tests)  is  denoted  for  this  work  as  the  structure  iden- 
tification procedure.  The  process  of  determining  the  coefficients  of  the  "identifica- 
tion" model  is  denoted  as  the  parameter  identification  procedure. 

2.5  IDENTIFICATION  OF  PARAMETERS 


The  extraction  of  aerodynamic  derivatives  from  flight  data  has  received 
considerable  attention  during  the  last  three  decades  and  the  most  recent  efforts 
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are  given  in  refs.  [14-15] . Earlier  techniques  were  mostly  manual  or  analog 
requiring  subjective  judgment  by  operators  [16] . These  methods  are  suitable 
for  simple  linear  systems  under  very  ideal  conditions.  A large  amount  of  data 
could  not  be  processed  because  of  lack  of  techniques  which  would  work  for  nonlinear 
system  models  in  the  presence  of  process  and  measurement  noise. 

More  recently,  with  the  availability  of  fast  digital  computing  machines 
and  efficient  computation  algorithms  many  powerful  digital  methods  have  been 
developed.  These  methods  are  usually  quicker,  do  not  require  subjective  judgment 
of  operators,  and  work  under  a variety  of  circumstances.  Examples  of  digital 
methods  are  various  equation  output  error  methods  [17-18] , output  error  methods 
[19-31]  , Kalman  filter/  smoother  approach  [10]  and  maximum  likelihood  technique 
[9.  32-34]. 

Output  error  methods , which  include  gradient  methods  and  modified  Newton- 
Raphson,  were  motivated  by  earlier  curve  fitting  techniques.  These  methods 
may  not  work  well  if  there  is  high  process  noise  or  if  the  weighting  matrix  in 
the  criterion  function  is  improperly  chosen.  The  equation  error  methods  minimize 
the  difference  between  the  left  hand  side  and  the  right  hand  side  of  the  state 
equations.  Examples  of  this  method  are  various  forms  of  least  square,  correlation 
methods,  and  instrumental  variable  approach.  Since  these  methods  do  not  account 
for  measurement  noise,  the  estimates  are  biased  in  the  presence  of  this  type 
of  error.  In  the  Kalman  filter /smoother  approach,  the  parameters  are  converted 
into  state  variables.  The  Kalman  filter  and  smoother  are  developed  for  this  new 
state  vector  starting  from  a priori  values  of  states  and  parameters  and  their  covari- 
ances. This  method  gives  biased  estimates  even  for  linear  systems  and  requires 
a priori  knowledge  of  measurement  and  process  noise  covariances. 

The  maximum  likelihood  technique  solves  some  of  the  problems  mentioned 
above.  By  considering  unknown  elements  of  process  noise  and  measurement  noise 
covariances  and  other  instrumentation  errors  as  parameters,  it  considers  these 
noise  sources  and  also  estimates  them,  if  not  known  a priori.  The  method  deter- 
mines parameter  values  which  maximize  the  likelihood  function  of  the  parameters, 
given  the  measurements  and  any  a priori  information.  The  likelihood  function 
has  the  same  form  as  the  conditional  probability  of  the  observations,  given  the 
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parameters.  The  maximum  likelihood  method  maximizes  the  probability  of  the 
observations.  It  is  customary  to  work  with  the  logarithm  of  the  likelihood  function. 
The  method,  as  applied  to  parameter  identification  in  nonlinear  dynamic  systems 
with  measurement  and  process  noise,  is  a combination  of  two  steps: 

1.  A Kalman  filter  for  estimates  of  the  state  and  its  covariance. 

2.  A Gauss-Newton  method  for  parameter  estimates,  and  associated 
covariances;  also  unknown  noise  statistics. 


The  details  of  the  above  steps  are  given  in  Appendix  A.  In  addition  to  estimating 
the  parameters  in  the  state  and  measurement  equations,  the  maximum  likelihood 
method  also  determines  the  covariance  of  errors  in  parameter  estimates . If  the 
model  whose  parameters  are  identified  is  a true  representation  of  the  system  in  the 
region  of  operation  and  the  sampling  rate  is  high,  the  maximum  likelihood  method 
gives  unbiased  estimates  of  parameters  for  long  data  records.  With  increasing 
amounts  of  data,  the  estimates  converge  to  their  true  value  almost  certainly.  It  can 
be  shown  that  the  technique  extracts  all  information  about  the  parameters  from 
data;  in  other  words,  the  method  is  efficient. 


2.6  INPUT  DESIGN 


Given  the  structure  and  parameters  of  a system,  it  is  possible  to  design 
experiments  which  enhance  the  accuracy  of  the  parameter  estimates.  Such  enhance- 
ment capability  is  extremely  important  for  design  of  flight  tests  either  from  pre- 
flight simulators  or  from  previous  flight  tests.  Analytically,  the  problem  of 
such  input  design  is  one  of  the  most  difficult  of  identification  technology. 


The  basic  idea  of  input  design  is  to  find  control  sequences  which  provide 
the  highest  possible  parameter  estimate  accuracy  in  the  shortest  time.  Recent 
advances  in  the  design  of  such  sequences  for  linear  systems  are  summarized 
in  Appendix  C of  this  report.  The  advances  achieved  in  linear  systems  are  much 
more  extensive  than  the  advances  for  nonlinear  systems.  Nevertheless,  the 
linear  techniques  may  be  selectively  applied  to  nonlinear  systems  for  various 
objectives.  (These  applications  are  discussed  in  Chapter  V.) 
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2 . 7 SUMMARY 


Procedures  for  validating  the  structure  and  parameter  estimates  are  denoted 
as  the  model  verification  process.  When  combined  with  instrumentation  and  input 
design  techniques,  the  integrated  parameter  identification  method  is  specified. 
Application  of  this  procedure  to  specific  problems  will  require  emphasis  of  some 
of  the  method's  specific  steps  and  placing  less  emphasis  on  others.  In  general, 
however,  all  steps  must  be  included  in  either  a qualitative  or  quantitative 
manner . 


Ill . MODEL  STRUCTURE  DETERMINATION 
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This  chapter  discusses  the  requirement  for  a method  of  estimating  which 
significant  aerodynamic  effects  cause  high  angle-of- attack  responses  (Section 
3.1).  This  requirement  is  usually  satisfied  by  basic  physical  considerations 
(e.g.  , the  equations  of  motion)  and  previous  wind  tunnel  and  flight  tests.  For 
the  complexity  of  the  high  angle-of-attack  regime,  however,  a subset  regression 
technique  is  presented  to  determine  an  aerodynamic  model  from  the  observed 
data  (Section  3.2)  . Examples  of  the  application  of  this  technique  show  that 
linear  models  are  frequently  useful  at  high  angle-of-attack  and  that  there  is 
useful  criterion  in  selecting  the  required  number  of  parameters  when  the  model 
is  nonlinear  (Section  3.3). 

3.1  FORMULATION  OF  AIRCRAFT  MODEL 

The  primary  requisite  of  any  parameter  identification  method  is  a model 

whose  parameters  are  to  be  estimated.  Preliminary  considerations  such  as  axis 

system,  equations  to  be  included,  and  type  of  aerodynamic  force  and  moment 

representations  are  discussed  in  Chapter  II.  More  detailed  considerations  include 

specifying  the  order  of  polynomials  to  represent  the  aerodynamics,  the  necessity 

for  including  terms  such  as  C and  C , and  the  determination  of  parameters  which 

ma  np 

should  be  known  a priori  to  achieve  the  best  identification  results.  The  latter 
aspect  is  based  on  the  fact  that  some  parameters  can  only  be  identified  by  extreme 
inputs  which  are  unsatisfactory  to  man  and  machine. 

3.1.1  Function  Representation 

The  main  problem  of  high  angle-of-attack  flight  data  reduction  for  stability 
and  control  coefficients  is  the  reconstruction  of  nonlinear  force  and  moment  equations. 
In  general,  this  may  require  estimation  of  more  parameters  than  is  computationally 
feasible.  In  addition,  the  requirement  to  identify  an  arbitrarily  large  set  of 
parameters  poses  significant  problems  of  identifiability . 

In  order  to  reduce  the  parameter  set  to  be  identified,  several  techniques 
may  be  applied.  Basic  to  all  of  these  techniques  is  the  requirement  for  "engineer- 
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ing  judgment"  to  postulate  a coefficient  from  and  evaluate  the  results  of  that 
assumption.  The  basic  techniques  are  outlined  as  follows: 

1.  Linear:  This  approach  seeks  to  identify  a set  of  linear  aerodynamic 

derivatives  for  many  flight  conditions.  The  objective  is  to  characterize 
a nonlinear  force  or  moment  variation  by  identification  of  function 
slopes  at  the  selected  flight  conditions . This  method  has  the  advantage 
of  associating  the  results  with  classical  stability  derivatives.  Converse- 
ly, it  does  require  tests  at  many  points,  some  of  which,  particularly 
in  the  high  angle- of- attack  regime,  may  not  be  susceptible  to  establish- 
ment of  a suitable  reference  state.  Difficulties  in  achieving  a predeter- 
mined steady  state  include  problems  of  controllability  and  buffet. 


2.  Quadratic  Function:  This  method  is  the  first  degree  of  complication 
beyond  the  linear  approach.  It  is  an  assumption  based  on  the  locally 
quadratic  behavior  of  airfoil  lift  characteristics  in  the  stall  regime. 

This  representation  has  been  used  for  at  least  two  decades  [35]  and 
is  still  widely  adopted  [36,37]  . The  representation  is  certainly  better 
in  nonlinear  regimes  than  the  linear,  but,  in  general,  will  have  greater 
error  at  the  limits  of  the  independent  variable  for  the  function  than 
higher  order  representations.  This  results  in  a problem  of  matching 
polynomial  segments  to  "build  up"  the  entire  polynomial. 


3.  Entire  Polynomial:  This  method  is  a further  extension  of  the  quad- 

ratic representation  which  allows  high-order  polynomials  to  represent 
the  desired  functions.  The  freedom  afforded  by  this  representation 
is  offset  by  the  difficulty  of  identifi ability  due  to  over-parameteriza- 
tion. In  addition,  this  representation  poses  strict  instrumentation 
accuracy  requirements  to  yield  satisfactory  reconstruction  of  high 
powers  of  the  function  independent  variable. 

For  many  modern  aircraft,  there  is  a comprehensive  series  of  wind  tunnel 
tests  which  provide  the  basis  for  design  evaluation  and  control  system  design. 

The  data  from  these  tests  may  also  be  used  to  formulate  functional  forms  which 
may  be  used  as  the  model  for  identification  from  flight  data.  Because  of  scaling 
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and  Mach  effects,  post  wind  tunnel  design  changes,  control  system  modifications, 
and  other  similar  effects  possibly  not  tunnel  tested,  these  functions  may  be  in 
error. 

The  most  desirable  goal  is  to  use  wind  tunnel  test  results  (and  any  other 
information  about  the  physics  of  the  maneuvers),  to  define  all  the  possibilities  of 
the  functions  for  the  forces  and  moments.  Then,  the  actual  data  obtained  would 
be  used  to  specify  which  of  these  functions  are  most  probable.  Identification  tech- 
niques would  then  be  applied  to  these  latter  functions . The  goal  is  to  identify  only 
those  polynomial  coefficients  which  are  required  to  reproduce  the  actual  force  or 
moment  characteristics.  To  this  end,  a fundamental  step  is  implemented  for  the 
integrated  parameter  identification  process. 

3.1.2  Model  Building  From  Data 

The  addition  to  the  integrated  parameter  identification  process  which  has 
been  implemented  for  high  angle-of-attack  application  is  diagrammed  in  Figure 
3.1. 

Flight  tests  are  conducted  for  particular  conditions  and  inputs.  The  data 
from  these  tests  is  then  used  as  the  basis  for  selecting  the  functional  forms  re- 
quii  ed  to  reproduce  the  data . The  data  is  passed  through  a filter  (denoted  as  the 
optimal  subset  regression  program).  This  filter  contains  suggested  variables 
which  might  be  required  to  reproduce  the  force  and  moment  characteristics 
which  generate  the  data.  These  variables  correspond  to  the  coefficients  of  powers 
of  the  independent  variables  (say,  a or  P)  . In  general,  the  number  of  terms 
allowed  is  dictated  by  expected  significance  from  a priori  considerations.  Of 
these  possible  required  terms,  the  filter  selects  the  most  significant.  The  model 
defined  by  these  significant  variables  is  then  passed  to  the  identification  program 
(Chapter  IV)  . 

The  polynomial  representations  chosen  by  this  approach  are  not  unique. 

Over  a finite  data  length  and  range  of  independent  variables , several  possible 
polynomials  may  fit  the  functional  form  adequately . For  example,  it  may  be  that 
the  two  forms 
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are  "equally"  good  representations  of  the  variation  of  pitching  moment  to  explain 
the  measurements  resulting  from  a certain  input.  The  optimal  subset  regression 
program  will  not  necessarily  favor  form  (a)  or  (b).  The  final  form  may  depend 
on  how  the  parameters  are  introduced  in  the  regression  equation.  The  coefficients 
in  form  (a)  and  (b)  are  different  and  are  chosen  such  that  both  functional  forms 
are  good  approximations  to  the  actual  variation  of  pitching  moment  with  angle- 
of-attack.  This  point  will  be  discussed  further  in  this  and  later  chapters  (e.g.  , 

IV  and  V)  . 


In  order  to  further  detail  the  use  and  limitations  of  this  program,  it  is 
necessary  to  discuss  the  basis  of  the  algorithm  formulation. 


3.2  SUBSET  REGRESSION  TECHNIQUES 


The  most  fundamental  approach  to  parameter  identification  is  historically 
through  the  original  least  squares  theory  of  Gauss.  Extensions  of  this  basic 
theory  include  modern  optimal  control  and  filter  theory,  as  well  as  the  widely- 
known  Kalman  filter  (a  sequential  least  squares  filter)  . The  extension  discussed 
here  is  conceptually  different  from  the  Kalman  filter  algorithm. 


The  theory  of  subset  regression  is  usually  found  in  advanced  texts  dealing 
with  statistical  inference.  The  mathematical  formulation  of  the  theory  is  reviewed 
in  Appendix  B,  and,  of  necessity,  is  directed  only  to  the  application  used  for 
this  work.  The  application  of  this  algorithm  to  aircraft  parameter  identification 
is  not  known  to  have  been  previously  reported. 


3.2.1  Basic  Concept  of  Subset  Regression 


Basically,  the  Kalman  filter  is  a linear  sequential  least  square  estimator  which 
provides  an  updated  minimum  variance  estimate  of  the  state  of  a dynamic  system  based 
on  observations  of  the  system  output,  y.  The  recursive  nature  of  the  algorithm  makes 
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it  possible  to  incorporate  the  information  from  the  most  recent  observation  without 
repeating  the  computations  already  made. 


Subset  regression  is  also  a sequential  procedure  except  that,  instead  of 
using  new  data  points  to  reduce  estimation  error  (over  time),  it  iteratively  con- 
verges on  a set  of  parameters  which  reduce  the  least  square  error  over  the  entire 
data  span.  The  resulting  set  of  parameters  constitutes  an  optimal  subset  of  all 
possible  parameters  which  significantly  correlate  the  input  and  output  data. 

The  principal  aspects  of  the  subset  regression  procedure  are  based  on 
the  decomposition  of  the  observations 


y = X0  + e 


X101  + X202  + E 


(3.1) 


(3.2) 


where  0 is  the  complete  set  of  parameters  and  0^  and  ©£  are  subsets  of  0.  The 
first  step  of  the  procedure  is  to  estimate  0^,  ignoring  0^, 
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ei=  (XJXj)  ^Xjy 


(3.3) 


The  estimated  residuals  are 


y = y - y 


(3.4) 
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where  ,y//  = [I  - X^(X^X^)  X^  ] . The  second  step  is  to  treat  y as  a new  observation 


and  obtain  the  least  square  estimate  for  ©_, , 


62  = (XTX2)-Ixjy 
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It  may  be  shown  [38]  that  this  technique  not  only  accounts  for  the  influence  of 
on  y,  but  also  on  y,  to  the  extent  that  and  X ^ are  correlated. 

The  implementation  of  the  algorithm  adds  and  deletes  variables  to  a particu- 
lar model  in  an  iterative  manner.  Estimates  of  previously  ignored  parameters 
are  incorporated  and  evaluated  by  two  criteria,  as  follows: 

1.  Of  all  possible  variables  0,  is  0.  the  most  highly  correlated  with  y 
of  variables  not  in  the  regression? 

2.  If  0.  is  added  to  the  regression,  is  its  contribution  to  the  "fit"  signifi- 
cant relative  to  variables  ••  • which  have  already  been 

used?  Does  the  significance  of  0.  , ....  ,0.  diminish  because  0.  is 
& l-l  l i 

included? 

These  questions  are  answered  within  the  framework  of  statistical  hypothesis 
testing  (discussed  in  Appendix  B)  . A generalized  flow  chart  is  shown  in 
Figure  3.2  . Starting  with  a list  of  possible  variables,  the  algorithm  enters  the 
variable  with  the  highest  correlation  (partial)  to  the  observations  y.  The  contri- 
bution of  this  variable  to  reducing  fit  error  is  made,  (5)  , (2)  , and  a new  variable 
entered,  (3),  Figure  3.2.  Subsequent  tests  add  and  delete  variables  to  improve 
the  "fit" . The  final  subset  of  0 which  results  bom  the  procedure  is  one  which 
falls  within  confidence  bounds  set  by  the  user  (for  example,  95%  or  99%). 

In  general,  the  subset,  0^,  of  0 is  not  unique.  The  particular  parameter  of 
the  final  subset  will  depend  on  the  order  in  which  they  are  listed.  This  occurs 
because  the  sum  of  the  squares  of  the  partial  correlations  of  each  variable  is  not 
equal  to  the  square  of  the  total  correlation  of  all  variables  on  y (ref.  Appendix  B)  . 
This  problem  m ay  be  reduced  by  using  orthogonal  variables  (which  could  meet 
the  sum  of  squares  equivalence)  or  by  forcing  in  variables  which  are  known  to  be 
required.  In  either  case,  a priori  correlations  must  be  introduced  which  depend 
on  knowledge  of  the  functions  sought.  The  implementation  used  for  the  purposes 
reported  here  includes  ability  to  force  some  variables  into  the  regression,  but  not 
reformulation  in  terms  of  orthogonal  variables. 
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Figure  3.2  Generalized  Flow  Chart  of  the  Subset  Regression  Algorithm 


The  algorithm  not  only  identifies  the  most  significant  parameters,  but  also 
finds  least  square  estimates  of  their  true  values.  In  general,  these  estimates  will 
be  in  error  (biased  due  to  measurement  noise  and  high  order  modeling  errors)  . 

As  such,  they  may  be  used  for  startup  values  of  the  maximum  likelihood  algorithm 
to  reduce  computation  time  and  improve  convergence,  but  must  not  be  considered 
as  final  estimates  themselves. 
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As  an  added  check  on  the  validity  of  the  determined  model,  the  program  com- 
putes the  residuals,  y-y,  on  the  final  pass.  This  provides  an  evaluation  of  the 
adequacy  of  the  model  (whose  values  yield  y)  compared  to  the  "true"  process  (whose 
values  give  y)  . Ideally,  these  residuals  will  be  white  Gaussian. 

3.2.2  Application  to  Nonlinear  Aircraft  Modeling 

The  parameters  with  which  the  high  angle -of- attack  applications  are  most 
concerned  are  the  coefficients  of  polynomial  expansions  of  the  aerodynamic  force 
and  moments.  The  reason  for  this  type  of  expansion  is  discussed  in  Chapter  II 
of  this  report. 

3.2.2. 1 Assumptions  on  Regression  Model 

The  following  assumptions  are  made  for  application  of  the  subset  regression 
algorithm: 

a.  Mass  and  inertias  of  the  aircraft  are  known.  This  assumption  is 
based  on  the  extensive  weight  and  balance  data  available  for  modern 
aircraft.  In  actual  test  conditions,  however,  there  may  be  marked 
differences  between  the  current  aircraft  configuration  and  the  base- 
line calculations  (e.g.,  unsymmetrical  loadings,  addition  of 
weapons  and  instrumentation).  This  assumption,  therefore, 
requires  evaluation  of  such  effects  to  determine  their  importance. 

b . '  Aircraft  rate  (rotary)  functions  are  assumed  to  vary  only  with 

ct  and  P . This  assumption  limits  the  model  structure  to  lower  Mach 

numbers  and  aircraft  loadings.  Specifically,  the  terms  C , C , 

md  n(3 

and  Cp>  are  not  considered  separately  from  C , C , and  C„  . 

^ n m n y. 

P q r r 

At  extreme  conditions,  this  assumption  is  questionable.  The  basic 
reason  for  allowing  use  of  this  assumption  is  that  the  data  base  [1] 
docs  not  have  the  aircraft  rate  and  velocity  vector  rate  terms  ex- 
plicitly separated.  In  general,  extreme  inputs  would  be  required 
to  separate  these  contributions  in  actual  flight  tests.  It  must  be 
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noted  that  the  regression  model  structure  can  be  easily  modified 
to  include  terms  such  as  C # if  required. 

P 

Aerodynamic  coupling  effects  due  to  combined  angle- of- attack 
and  sideslip  angles  are  limited  to  second  order  terms . This 
assumption  limits  the  cross -coup ling  terms  of  the  aerodynamic 
expansions  . The  expansions  in  terms  of  a alone  or  p alone  are 
unlimited,  however. 


d.  Measurements  of  aircraft  linear  and  angular  accelerations  are 

available.  This  assumption  is  the  result  of  the  regression  procedure. 
In  particular,  the  method  requires  measurements  of  p , q,  and  r, 
in  addition  to  the  linear  accelerometer  readings.  If  not  available, 
equivalent  combinations  of  roll,  pitch  and  yaw  rates  may  be  used. 

For  the  examples  considered  in  this  chapter,  the  measurements  are 
assumed  noiseless . This  assumption  is  dropped  in  Chapter  IV . 


3 . 2 . 2 . 2 Specific  Expansions 

The  general  expansion  for  any  specific  force  or  moment  coefficient  is 


C = Co(a,p)  +2C  (i)a(l)  + ZC 


p(])P  i j a«p0> 


where  a and  P are  reference  an gle-of- attack  and  sideslip  angle,  respectively.  The 
specific  coefficients  for  which  such  an  expansion  is  used  are  listed  in  Table  3.1. 


Primary  emphasis  is  placed  on  the  static  coefficient  structure  for  the 
modeling  tasks  here.  All  static  terms  could  be  minimally  expressed  as  linear 
combinations  of  (C.a  + C.p)n , where  n=l,2,3.  In  addition,  the  coefficients  of  Cm 
could  be  expanded  to  a9  and  C to  p5.  Control  effectiveness  coefficients  and 
dynamic  coefficients  were  generally  limited  to  first  order  expansions  in  a only. 
These  expansions  were  decided  upon  after  initial  experimentation  with  the  pro- 
grams . It  was  found  that  higher  than  first-order  terms  in  a were  usually  not 
required  on  the  basis  of  the  regression  analysis.  In  general,  these  control  effect- 
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TABLE  3.1 


I 


I i 

I 

p.  ' 

K i 


1 


] 


Coefficients  Which  Are  Expanded  for  High  Angle-of- Attack 
Aerodynamic  Models 


EQUATION 

STATIC 

COEFFICIENT 

CONTROL 

EFFECTIVENESS 

DYNAMIC 

COEFFICIENT 

u 

C (a,p) 

C 

C 

x6 

s 

X 

q 

V 

C (a,F) 

C C 

c 

a r 

yp 

w 

C (a,p) 

C 

c 

z 

z6 

s 

Z 

q 

p 

C£  (a,p) 

Cp  Cp 

6 6 
a r 

C£  C£ 

P r 

q 

C (a,p) 

C 

C 

m 

m6 

a 

m 

q 

r 

C (a,p) 

C C 

c c 

n 

n6  n8 
a r 

n n 

P r 

iveness  and  dynamic  terms  are  difficult  to  identify  from  acceptable  inputs,  and 
first-order  a expansions  were  considered  acceptable  . 


An  important  aspect  of  the  regression  program  is  that  the  coefficients  are 
not  serial.  Thus,  if  intermediate  coefficients  were  not  necessary  (e.g.  , such  as 
Ca2p),  they  are  not  included.  (The  maximum  likelihood  program  detailed  in 
Chapter  IV  will  automatically  eliminate  terms  not  specified  by  the  regression.) 


3. 2. 2. 3 Equation  Error  Formulation  of  Dynamic  Responses 

The  basic  objective  of  the  regression  program  is  to  determine  the  aero- 
dynamic coefficients.  Implicit  is  the  requirement  that  all  other  variables  (in 
particular,  the  aircraft  states)  be  known.  The  general  equations  of  motion 
(Appendix  E)  may  be  written  in  the  form 

x = f(x,u,0)  + e (3.10) 


__  , ......  ...  . , 


: 
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The  subset  regression  is  applied  to  Eqs  . (3.1)  equation-by-equation. 

The  left  hand  sides  of  Eqs . (3.1)  are  shown  in  Table  3.2.  Note  that  the  dynamic 

pressure  is  included  in  this  left  hand  side.  The  variables  ILU,  ELV , ELW, 

WLP,  ELQ,  ELR  are  calculated  from  the  measured  values  of  the  aircraft  state, 

as  detailed  in  Table  3.2.  These  "pseudo-measurements"  are  formed  from 

the  actual  measurements  such  as  p^,  q^,  (the  measured  values  ofp,  q,  r) 

and  accelerometer  measurements  (x  . v . z ) 

m 7 m m 

For  this  application,  the  terms  ELU . . . ,ELR  are  actually  calculated  in  the 
simulation  program  used  to  generate  the  data.  In  actual  test  data  processing , a 
separate  program  would  have  to  be  used  to  formulate  these  terms.  It  is  also 
required  to  substitute  measured  values  of  a,p,  and  the  control  deflections  into  the 
program  to  calculate  the  observations  of  "independent"  variables  (h(a,p)  of  Eq. 

(3.13)).  This  latter  requirement  introduces  a modeling  error  into  the  estimation 
scheme . 

3.3  EXAMPLE  OF  SUBSET  REGRESSION  APPLICATION  TO  SIMULATED  DATA 

3-3,1  Lateral-Directional  Response  Application  of  the  Subset  Regression 

Method 

One  of  the  first  steps  in  evaluating  the  subset  regression  program  was  the 
lateral  directional  response  analysis  . This  case  was  considered  to  be  more 
important  and  difficult  than  the  longitudinal  tests . Most  of  the  reported  incidences 
of  degraded  stall/post-stall  responses  occur  in  roll  or  sideslip,  and  multiple 
non  linearities  characterize  the  regime. 

A number  of  small  perturbation  tests  were  performed  for  this  case.  It  was 
determined  that  severe  or  complicated  control  time  histories  were  not  warranted 
for  flight  tests  at  high  angle-of- attack  (due  to  the  basically  unstable  character  of 
responses  in  the  regime)  . The  fundamental  aspect  of  these  examples  was  to 
determine  the  adequacy  of  a linear  model  at  several  angles-of-attack , as  opposed 
to  a polynomial  model  requiring  more  complicated  inputs.  These  results  are  shown 
in  Tables  3.3  and  3.4. 


37 


TABLE  3.3 


Aileron  Doublet  (+5°  over  2 second  period,  10  sec  data  length) 


ANGLE  OF 
ATTACK 

MODEL  REQUIRED  TO  EXPLAIN  AT  LEAST  99%  OF 
Cg  VARIATION 

5° 

C£  = C£  P + C£  6a  + Cfi  p 
P 6a  p 

15° 

C£  = C£  P +C£  6a  + C£Z.2a2P2 
B *6  Vb 

a r 

25° 

C£  = C£  P + Cfi„  6a  + C£  P + C£  r 

P 8a  p r 

TABLE  3.4 

Rudder  Doublet  (+5°  over  2 second  period,  10  second  data  length) 


ANGLE  OF 
ATTACK 

MODEL  REQUIRED  TO  EXPLAIN  AT  LEAST  99%  OF 

C VARIATION 
n 

5° 

C=C  6+C  p+C  r 

n nR  r n„  n 

6r  P r 

15° 

C = C 8+C  p+C  r + C p 

n n-  r n„  n n ^ 

6r  P r p 

25° 

C=C  6+C  p+C  r 

n nR  r n„'  n 

o P r 

r ' 

Note  that  there  is  roll-pitch  coupling  at  15°  angle-of-attack  for  the  aileron 

doublet  example.  This  departure  from  linearity  occurs  because  of  the  pitch 

instability  at  that  point.  In  each  case,  Cg  =0,  which  indicated  zero  rolling  moment 

o 

at  trim.  In  addition,  the  wind  tunnel  results  indicated  Cg  - 0 at  a = 15°  which  was 

p 

substantiated  by  model  structure  determination  program.  At  least  98%  of  the 

moment  variation  in  each  <xme  was  explained  by  using  Cp  , Cp  and  C0  , with  the 

P 6 o 

ci 

other  terms  being  much  less  significant.  Again,  in  each  case  C =0,  indicating 

o 

a zero  yawing  moment  at  trim.  The  C coefficient,  as  verified  by  the  known  data 

P 

changed  sign  (from  positive  to  negative)  between  a = 15°  and  a = 25°.  As  in  the 
aileron  doublet  case,  at  least  98%  of  the  yawing  moment  variation  could  be  explained 
by  using  simply  and  . A linear  model  seems  quite  adequate  in  this 

6r  P r 

case  where  the  critical  sign  change  in  C was  found  by  this  approach. 

np 

The  ability  of  a linear  model  to  explain  perturbation  input  aircraft  re- 
sponses, even  in  a "nonlinear"  regime,  has  been  discussed  in  Section  3.1.1. 

The  linear  mode]  is  here  shown  to  indeed  be  adequate  as  a point-wise  definition 
of  the  aerodynamic  model.  This  indicates  the  possibility  of  using  linear  input 
designs  for  high  angle-of-attack  regimes,  a subject  elaborated  on  in  Chapter  V. 


This  conclusion,  however,  must  be  approached  with  caution.  At  high 
angles  of-attack,  it  may  be  difficult  to  hold  trim  and  to  conduct  even  perturbation 
maneuvers  at  constant  angle-of-attack.  Not  only  inherent  instability,  but  also 
buffet,  may  make  such  inputs  difficult  at  best.  The  important  conclusion  which 
can  be  attained  is  that  the  regression  program  can  be  used  to  delineate  those 
regions  where  linear  or  nonlinear  aerodynamic  modeling  is  required. 


3.3.2  Post-Regression  Subset  Specification 

The  subset  regression  approach  is  an  elegant  and  efficient  means  of  model 
structure  estimation.  As  with  most  data  processing  aids,  however,  its  use  must 
be  moderated  with  engineering  judgment.  Extensive  use  of  the  technique,  as  with 
any  analytical  approach  repeatedly  applied  to  physical  problems,  allows  formula- 


tion  of  certain  guidelines  which  facilitate  such  judgment.  Such  guidelines  usually 
have  a basis  in  the  theory,  but  are  not  easily  quantified. 


One  such  guideline  is  the  Akaike  final  value  prediction  theorem  [39].  This 

criterion  states  that,  for  certain  types  of  systems,  there  is  an  optimal  number  of 

parameters  which  describe  a model  structure.  A plot  of  this  criterion  function 

versus  the  number  of  parameters  has  a unique  minimum  (i.e.,  in  a quadratic 

sense)  which  is  at  the  optimal  number  of  parameters  [9] . It  is  also  shown  that 

fit  error  (the  residuals  of  the  estimated  versus  actual  output  time  history)  is  not 

2 

satisfactory,  since  it  (as  also  measured  by  R , the  multiple  correlation  coefficient 
of  Appendix  B) , as  a criterion,  does  not  possess  a unique  minimum. 


Fit  error  criteria  tend  to  approach  an  asymptotic  value  as  the  number  of 
parameters  is  increased.  The  eventual  insensitivity  of  the  fit  error  to  an  increase 
in  the  number  of  parameters  (as  parameters  are  added)  , is  due  to  the  continuing 
reduction  of  degrees-of -freedom . * Fundamentally,  when  the  number  of  data  points 
equals  the  number  of  parameters,  the  regression  curve  passes  exactly  through 
these  points.  No  further  improvement  is  then  possible.  Noise  will,  of  course, 
allow  more  parameters  to  be  entered  (since  point-by-point  fit  no  longer  occurs)  , 
but  these  parameters  tend  to  fit  the  noise,  not  the  process. 

The  Akaike  criterion  is  1 function  of  fit  error,  but  also  weights  this  error 
with  the  number  of  degrees  of  freedom.  For  the  regression,  a similar  function  of 
fit  error  is  the  F-ratio,  the  ratio  of  "fit  goodness"  to  fit  error,  weighted  by  the 

degrees-of-freedom . The  concept  of  F-ratio  is  detailed  in  Appendix  B,  and  its  im- 
portance to  this  work  is  significant. 


That  the  F ratio  is  a measure  of  the  "optimal  number"  of  parameters  is  shown 

in  Figure  3.3  , corresponding  to  an  example  to  be  discussed  in  more  detail  in 

Chapter  VI.  Its  placement  here  is  to  emphasize  the  need,  based  on  a criterion 

2 

other  than  fit  error,  to  evaluate  model  suitability.  As  shown,  R is  a inonotoni- 
cally  increasing  function  of  the  number  of  parameters.  The  F-ratio,  however, 
has  a maximum  with  seven  parameters  (obtained  after  deleting  one  variable  of  small 

* 

Number  of  observations  less  number  of  parameters. 
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significance)  . 

In  general,  it  has  been  found  that  the  subset  F-ratio  will  have  local  maxima, 
beyond  the  first,  as  parameters  are  added.  This  result  is  shown  in  Figure  3.4  . 
The  first  general  criterion  used  for  selecting  parameters  to  be  identified  is  to 
delete  all  parameters  included  in  the  regression  past  the  first  maximum. 

A second  guideline  to  further  optimize  the  regression  subset  for  maximum 
likelihood  application  is  to  select  only  parameters  whose  contribution  is  a certain 
percentage  of  the  most  significant  parameter.  This  is  discussed  in  Chapter  VI. 

3.4  SUMMARY  OF  CAPABILITY  OF  THE  SUBSET  REGRESSION  PROCEDURE 

This  chapter  has  presented  the  basic  idea,  implementation  details,  and 
example  results  of  the  model  determination  subset  regression  method.  Extensive 
use  of  this  program  will  be  discussed  in  more  detail  in  Chapter  VI.  At  this  point, 
a review  of  the  program  features  will  be  emphasized. 

1.  The  program  gives  a set  of  parameters  which  approximately  describe 
the  input-output  relationship.  The  set  of  parameters  may  not  be 
unique , however  . 

2.  The  procedure  produced  biased  parameter  estimates  in  the  presence 
of  measurement  noise.  The  method  is  based  on  a least  squares  prin- 
ciple, and  measurement  errors  are  well-known  to  produce  biased 
estimates  for  this  technique. 

3.  The  method  can  be  used  for  input  design  for  nonlinear  as  well  as 
linear  systems. 

4.  The  method  is  extremely  rapid  for  producing  a model  structure 
estimate.  This  has  been  found  from  extensive  computation  on  a 
UN1VAC  1108. 
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5. 


The  method  can  be  used  to  provide  a priori  (though  biased)  estimates 
of  parameters . This  follows  since  the  program  can  basically  be 
interpreted  as  an  equation  error  start-up  procedure. 


The  robustness  of  the  program  to  process  and  measurement  noise,  and  its 
operational  use,  are  discussed  in  the  following  chapter. 

This  model  determination  technique  is  the  first  step  of  the  method  for  iden- 
tifying parameters  from  flight  test  data  in  the  high  angle- of- attack  regime.  At 
this  point  in  the  discussion  of  the  method,  a procedure  has  been  established  for 
isolating  the  cause-effect  relationships  necessary  for  modeling  the  aircraft  re- 
sponse. Accurate  quantification  of  these  models  is  performed  by  the  second  stage 
of  the  method — the  maximum  likelihood  procedure. 


IV . MAXIMUM  LIKELIHOOD  METHOD  FOR  NONLINEAR  SYSTEMS 


An  overall  view  of  the  parameter  identification  process  for  aircraft  stability 
and  control  derivative  determination  was  presented  in  Chapter  II.  In  the  high 
angle -of- at  tack  flight  regime,  where  the  aircraft  motions  are  nonlinear,  the  first 
step  in  the  identification  procedure  is  the  determination  of  the  model  structure 
which  adequately  describes  the  relationship  between  control  inputs  and  measured 
aircraft  response.  The  subset  regression  method  of  Chapter  III  identifies  the 
model  and  gives  a first  pass  at  the  values  of  unknown  parameters  in  the  model. 

This  is  computationally  simple  and  quick.  The  final  step  is  the  use  of  the  maximum 
likelihood  method  to  refine  the  model  structure  and  to  improve  parameter  estimates 
so  that  they  are  unbiased  and  have  minimum  variance,  even  in  the  presence  of 
process  and  measurement  noise. 

In  this  chapter,  the  concept  of  the  maximum  likelihood  method  is  reviewed 
(Section  4 . 1)  . The  specific  maximum  likelihood  features  used  for  this  project 
are  detailed  (Section  4.2),  and  related  to  the  detailed  mathematical  equations 
upon  which  these  features  rest  (Appendix  A)  . Computational  considerations 
are  presented  (Section  4.3)  and  an  example  output  is  given  for  illustration  (Section 
4.4)  . 

4 . 1 INTRODUCTION 


The  maximum  likelihood  method  is  described  in  Appendix  A and  is  illus- 
trated in  Figure  4.1.  Conceptually,  this  technique  can  be  summarized  as  follows: 

"Find  the  probability  density  functions  of  the  observations  for  all 
possible  combinations  of  unknown  parameter  values.  Select  the  den- 
sity function  whose  value  is  highest  among  all  density  functions  at 
the  measured  values  of  the  observations.  The  corresponding  para- 
meter values  are  the  maximum  likelihood  estimates." 
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FLIGHT  TEST  DATA,  WIND  TUNNEL  VALUES  OF 


AERODYNAMIC  PARAMETERS 


Figure  4.1  Flowchart  of  Maximum  Likelihood  Identification  Program 
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Suppose  9 can  take  three  possible  values:  0-^,  G^,  and  G^.  Let  the  probability 
density  functions  of  observations  z for  these  three  values  of  9 be  as  shown  in 
Figure  4.2.  Then,  if  the  actual  observation  is  z,  6^  is  the  maximum  likelihood 


estimate  of  8.  In  practice,  starting  from  a priori  estimates,  the  parameters  are 
updated  so  that  the  value  of  the  resulting  density  function  at  the  observations 
increases  monotonically . 

The  concept  implementation  is  shown  in  Figure  4.2.  An  extended  Kalman 
filter  is  used  with  fixed  parameter  values  to  determine  state  and  covariance  time 
histories.  The  likelihood  function,  which  has  the  same  form  as  the  density  function 
of  the  observations,  depends  on  innovations  and  their  covariance.  The  first  and 
second  gradient  of  the  likelihood  function  are  calculated  by  propagating  equations 
governing  the  sensitivities  of  state  and  covariance  called  sensitivity  equations. 
Gradient  procedures  are  used  to  determine  step  size  and  update  the  parameters  . 
This  is  repeated  until  convergence  occurs  . 

The  most  important  asset  of  the  maximum  likelihood  method  is  its  versatility. 
The  next  sections  describe  how  this  technique  is  used  for  estimating  parameters 
in  nonlinear  aircraft  models  at  high  angles-of-attack . 

4.2  MAXIMUM  LIKELIHOOD  TECHNIQUES 

The  maximum  likelihood  method  is  an  inherently  powerful  technique  for  the 
study  of  parameter  identification  in  the  high  angle-of-attack  flight  regime.  The 
aircraft  motions  in  this  region  are  highly  nonlinear  and  there  is  no  single  model 
to  adequately  describe  the  motions.  Here  the  method  is  simplified  so  that  more 
attention  is  given  to  the  identification  and  related  problems  of  identifiability 
of  parameters  in  the  model  structure  determined  by  techniques  of  Chapter  III. 
Certain  simplifications  which  are  made  disregard  effects  which  are  only  secondary 
to  the  above  study.  These  modifications  are  detailed  in  the  following  paragraphs. 


"observed 


Figure  4.2  Maximum  Likelihood  Estimates 


In  the  presence  of  process  noise,  parameter  estimates  may  be  biased  unless 
estimates  of  the  process  noise  statistics  themselves  are  made.  To  estimate  the 
process  noise  statistics,  it  is  necessary  not  only  to  propagate  the  state  and  sensi- 
tivity for  all  unknown  parameters,  but  also  to  propagate  the  covariance  of  the  state 
estimation  error  and  its  sensitivities  . This  is  an  increase  in  computational  require 
ments  over  the  case  where  process  noise  statistics  are  not  identified.  The  ability 
of  the  maximum  likelihood  algorithm  to  identify  process  noise  statistics  has  previ- 
ously been  demonstrated  [9]  . The  resolution  of  the  problems  of  ovei  parameter- 
ization is  the  central  objective  of  this  work.  Therefore,  the  identification  of  process 
noise  statistics  is  not  included  in  this  implementation  of  the  maximum  likelihood 
algorithm.  As  will  be  shown,  however,  parameter  estimates  are  accurately  ob- 
tained, even  with  process  noise. 

The  algorithm  equations  are  detailed  in  Appendix  A of  this  report.  This 
appendix  contains  not  only  the  fundamental  algorithm,  but  also  the  various  analyti- 
cal forms  which  the  algorithm  assumes  under  various  simplifications.  In  particular, 
the  algorithm  used  for  the  results  of  this  work  is  obtained  by  noting  that 


xOlt.i)  = xCt) 


(4.1) 


The  implementation  of  the  maximum  likelihood  method  used  here  is  quite 
sophisticated.  Some  of  the  main  features  of  the  computer  program  used  to  study 
parameter  identification  at  high  angles-of-attaclc  are; 

(a)  The  program  incorporates  a detailed  instrument  model.  The  measure- 
ments are  assumed  to  have  bias,  scale  factor  errors,  and  random 
noise.  For  example,  the  measured  value  of  pitcHrate  is 


q = (1+k  )q  + b + w 

Mm  q q q 


(4.2) 


The  scale  factor  k , bias  b and  the  spectral  density  of  w can  be 

q q q 

identified  or  given  preassigned  values. 


(b)  In  most  implementations  , the  measurement  noise  covariar^ce  matrix 
is  estimated  by  . 


a i rT' 

E v(i)  v (i) 
i=l 


(4.3) 


where  v(i)  are  innovations.  In  general,  for  finite  data  length,  this 
would  give  an  estimate  of  the  measux-ement  noise  covariance  matrix 
with  all  nonzero  elements  and  no  special  structure.  In  other  words, 

Eq.  (4.3)  assumes  that  there  are  m(m+l)/2  unknown  elements  in 
measurement  noise  statistics.  If  the  structure  of  R is  known,  as 
is  usually  the  case  in  aircraft  applications,  the  number  of  unknown 
elements  in  R is  fewer  than  m(m+l)/2.  For  example,  if  the  measurement 
noise  in  different  channels  is  independent,  only  diagonal  terms 
of  R need  to  be  estimated.  The  use  of  Eq.  (4.3)  will  lead  to  an  over- 
parameterized  model  and  incorrect  parameter  estimates.  The  statisti- 
cally correct  estimate  of  the  measurement  noise  covariance  matrix 
for  independent  measurement  noise  sources  is, 
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(4.4) 


a 1 ( N -p 

R=  diag  Nl.2  v 

(c)  In  the  high  angle-of-attack  flight  regime,  the  longitudinal  and  lateral 
motions,  in  addition  to  being  nonlinear,  are  coupled.  This  requires 
the  solution  of  nine  nonlinear  differential  equations  in  u , v , w , 

p,  q,  r,  cp,  G and  v|/.  In  many  instances,  one  or  more  equations 
do  not  have  any  unknown  parameters.  The  program  is  written  so 
that  some  state  equations  can  be  discarded.  Since  the  aircraft  is 
fairly  well  instrumented,  there  are  measurements  of  all  state  variables. 
The  measurements  of  those  states,  whose  governing  equations  have 
been  discarded,  are  used  in  other  equations  instead  of  the  value 
of  the  states  obtained  by  integrating  these  equations . This  results 
in  a considerable  saving  in  computation  time  at  a slight  loss  in  accur- 
acy. This  feature  can  be  used  to  identify  longitudinal  and  lateral 
derivatives  separately  even  in  regions  where  these  two  motions 
are  coupled . 

(d)  As  explained  in  the  last  chapter,  the  aerodynamic  derivatives  are 
expressed  as  polynomial  functions  of  angle-of-attack  and  sideslip 
angle.  No  single  expansion  with  a reasonable  number  of  terms 
describes  most  of  these  derivatives  adequately  in  the  entire  range 
under  consideration.  Therefore,  different  expansion  terms  are  re- 
quired in  diffei'ent  operating  regions.  The  computer  program  can 
handle  all  possible  expansions  of  aerodynamic  coefficients.  For 
example,  the  equation  governing  pitch  rate,  excluding  nonlinear 
dynamic  terms,  is, 


1 


] 


■ ' — ■ — MM  ■ 


, a dCm  (a,P)  dC  (a, p) 

dF  £ = ^ + Cm  (a'P)  ^ + ~ 67 

q 

dC  («,P) 

ill  p 

+ ^ 6 

dc  s 


(4.6) 


If  the  expansion  for  in  a particular  case  is 


k 2 

C (a,p)  = 2 SC..  aV 
m i=0  j=0  1J 


(4.7) 


where  some  C./s  are  zero,  then 


dc 


dC  (a,p)  k 2 ( 3C..  . . . . . . . ) 

= 2 2 +iC..  a1_1p^  — + jC..  a'pJ'V 

i=0  j=o  l 9c  ^ P 9c  J 13  P \ 


(4.8) 


and 


3C..  I = 1 if  c is  C. . 

LL  < l3 

9c  f = 0 otherwise 


(4.9) 


This  implementation  of  the  maximum  likelihood  method  can  compute 
the  gradient  of  C^Ca.p)  and  other  aerodynamic  derivatives  with 
respect  to  any  parameter  c for  chosen  k and  2 and  any  zero. 

(e)  The  computer  program  determines  the  covariances  of  estimates  of 
states,  measurements  and  aerodynamic  coefficients  as  a function 
of  time . Let  the  state  and  measurement  equations  be 


x = f(x,u  , 0,  t) 


(4.10) 


y = h(x,u,0,t)  + v 


(4.11) 
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The  error  in  state  estimate  is 


- _ 9x  ^ 

x " W 9 


The  covariance  of  state  estimation  error  is 


>,,-1  (dx)1 

ae  M [do) 


where  M is  the  information  matrix,  the  inverse  of  which  is  the  Cramer 
Rao  lower  bound  on  parameter  error  covariances. 

The  equations  for  computing  state  sensitivities  are  given  in  Appendix 
A.  The  output  of  the  model  is 

y = h(x,u,0,t) 


The  covariance  of  output  error  is  approximat 


-♦&K1  &)*♦» 


The  correlation  between  0,  x and  v is  neglected.  The  covariance 
of  error  in  aerodynamic  coefficients  is  computed  in  a similar  fashion. 


(f)  It  is  possible  to  fix  any  parameter  at  a preassigned  value  or  constrain 
it  to  stay  within  a region. 

The  expansion  terms , which  are  identified  to  be  significant  for  an  aero 
dynamic  coefficient  during  model  building,  do  not  completely  explain  its  variation 
with  a and  p.  Therefore,  the  model,  whose  parameters  are  identified,  lias  residual 
terms  called  modeling  errors.  These  modeling  errors  produce  biased  parameter 
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estimates.  The  bias  in  the  estimates  depends  on  the  amount  of  modeling  errors 
and  how  changes  in  parameters  can  explain  this  error.  Bias  also  results  from 
inaccurate  instrument  modeling  or  fixing  parameters  at  incorrect  values.  The 
expected  values  of  these  biases  do  not  depend  upon  the  amount  of  data.  On  the 
other  hand,  the  covariance  of  parameter  estimates  decreases  with  increasing 
data  size.  In  a particular  run,  it  is  worthwhile  to  make  simplifications  which 
produce  biases  in  parameters  much  smaller  than  the  corresponding  standard 
deviations.  This  could  lead  to  a considerable  saving  in  computation  time  for 
an  insignificant  loss  in  accuracy.  For  instance,  simplified  instrument  models 
may  be  adequate  for  short  flight  testing  time. 

Information  Aggregation 


Any  available  a priori  information  about  parameters  can  be  incorporated 
into  the  identification  program.  This  is  done  by  finding  the  weighted  mean  of 
the  a priori  parameter  values  and  the  output  of  the  identification  program,  the 
weights  varying  inversely  as  the  corresponding  covariances.  This  has  been  done 
before  T9]  . A similar  approach  can  be  used  to  combine  estimates  from  different 
identification  runs  and  thus  make  full  use  of  available  data.  These  runs  give 
estimates  of  aerodynamic  stability  and  control  derivatives  over  different  ranges 
of  a and  p,  together  with  the  corresponding  covariances.  Since  these  coefficients 
are  approximated  by  low  order  polynomials,  the  coefficients  of  the  polynomial 
provide  a "good"  fit  over  the  range  of  a and  p covered  by  the  experiment.  It 
is  not  possible  to  use  these  polynomials  for  extrapolation.  Therefore,  the  informa- 
tion about  the  value  of  an  aerodynamic  coefficient  obtained  from  all  experiments 
at  discrete  values  of  a and  p can  be  combined . Suppose  that  k different  experiments 

give  estimates  C.  with  variances  o.  for  the  value  of  C at  a . Then,  the  "best" 
i i mo 

estimate  of  C (a  ) would  be 
m o 


A 

c 

m 


(ao)  = 


k 

Z 

i=l 

k 

Z 

i=l 


C. 

1 

a 

i_ 

J_ 

a. 

i 


(4.15) 
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This  can  be  done  at  many  discrete  points  to  obtain  a continuous  curve. 
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4.3  COMPUTATIONAL  ASPECTS  OF  MAXIMUM  LIKELIHOOD  ESTIMATES 


! 


; 


i 


The  application  of  maximum  likelihood  identification  in  practice  requires 
accurate  computational  algorithms  to  maximize  the  likelihood  function.  It  is 
not  uncommon  to  find  situations  where  the  likelihood  surface  has  multiple  maxima, 
saddle-points,  discontinuities  and  a singular  Hessian  in  the  parameter  space.  The 
application  of  the  steepest  descent  method  leads  to  an  extremely  slow  convergence 
rate,  and  the  straightforward  application  of  the  Newton-Raphson  and  the  Gauss- 
Newton  methods  may  lead  either  to  no  convergence  or  convergence  to  wrong 
stationary  points.  From  a statistical  viewpoint,  only  the  absolute  maximum  of  the 
likelihood  function  provides  an  unbiased,  consistent,  and  efficient  estimate.  Thus, 
it  is  important  to  locate  the  true  maximum  of  the  likelihood  function.  With  the 
present  computation  methods  of  nonlinear  programming,  in  general,  this  could 
be  an  extremely  difficult  and  time-consuming  task.  In  aircraft  applications,  a priori 
information  about  the  parameters  is  used  to  choose  a good  starting  value  and 
impose  physically  meaningful  constraints  on  the  parameters.  This  increases 
the  chances  of  approaching  the  absolute  maxima. 

4.3.1  Identifi ability  Problems  and  Solutions 

Anomalies  also  occur  in  the  likelihood  function  due  to  inadequate  model 
specification  and  parameterization.  Some  of  the  anomalies  that  lead  to  numerical 
difficulties  are: 


Discontinuities  and  Singularities:  There  are  examples  where  the 

likelihood  function  increases  up  to  a maximum  value  and  then  falls 
immediately  to  minus  infinity.  In  other  cases,  the  likelihood  function 
may  rise  suddenly  to  infinity  for  some  particular  parameter  values. 
These  problems  can  often  be  avoided  by  either  changing  the  model 
or  changing  the  parameterization. 

Singular  or  Nearly  Singular  Information  Matrix:  This  is  one  of 

the  most  common  problems  in  parameter  identification  since  it  is 
difficult  to  determine  a priori  whether  a given  sample  has  adequate 
information  for  estimating  all  the  parameters  in  the  model.  Generally, 
the  tendency  in  practical  applications  is  to  include1  all  the  parameters 
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about  which  there  is  some  uncertainty.  This  leads  to  overparameteri- 
zation which  in  turn  produces  a singular  or  a nearly  singular  informa- 
tion matrix.  The  likelihood  surface  in  turn  has  long,  curved  and 
narrow  ridges  along  which  the  convergence  rate  is  extremely  slow. 
Under-parameterization  does  not  really  provide  the  solution  since 
it  may  lead  to  spurious  local  maxima  and  saddle-points. 


The  basic  iteration  in  gradient-type  nonlinear  programming  methods  is 


e . i = e.  - p.  r7  g. 

i+i  i i bi 


(4.16) 


where  0.  is  the  parameter  vector  at  the  i 1 iteration,  g.  is  a vector  of  gradients 

^ *th 

of  the  negative  log-likelihood  function  J(0)  , i.e. , the  j 1 component  of  g.  is  given 


_ v*  ) T„-l  9v  , T„-l  9R  ^-1 

8i(j)  ^ j 90.(j)'iV  R aG.(j)  R V 


, , ^ /n-l  3R  ' 

+ iTr\R  aoTTf) 


(4.17) 


R^  is  an  approximation  to  the  second  partial  matrix 


and  is  a scalar  step  size  parameter  chosen  to  ensure  that  J(0^+p  < J (0. ) - e 
where  e is  a positive  number  that  can  be  chosen  in  a variety  of  ways.  In  the 
Gauss-Newton  method  used  here,  the  matrix  of  second  parti als  is  approximated 
by 


Vi.x)  -mi0.w  wffc 


, i T n-l  9R 

tlTTR 


,-1  9R 


(4.18) 
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where  M is  the  approximation  of  R.  The  first  partials  of  innovations  v and  their 
covariance  can  be  computed  by  solving  linear  difference  equations  (i.e. , the 
sensitivity  equations).  Notice  tli at  in  this  approximation,  computation  of  second 
partials  of  state  and  covariance  is  not  required . 

There  are  several  problems  with  this  method.  As  long  as  M.1  is  positive 
semi- definite , it  is  possible  to  find  a step  which  will  reduce  the  cost.  As  mentioned 
before,  most  of  the  problems  arise  when  the  information  matrix  is  almost  singular 
(i.e. , the  eigenvalues  are  spread  far  apart)  , which  is  almost  always  the  case  in 
high  angle-of- attack  problems  . This  near-singularity  of  the  information  matrix 
manifests  itself  in  several  ways: 

(a)  During  the  computation  of  the  inverse  of  the  information  matrix, 
the  positive-definiteness  may  be  lost  because  of  round-off  errors. 

(b)  The  step-size  may  be  very  large  in  those  directions  which  correspond 
to  small  eigenvalues  of  the  information  matrix.  Let 


M = T k.v.v. 
i=l  1 1 1 


(4.19) 


M"1 - y 1 T 

M =2.  -r—  V.V. 

. , X.  11 
1=1  1 


(4.20) 


“ 1 T 

a.E,rviviS 

1=1  i 


m (viT  g) 

- E c<4_J  v. 

i=l  \ 1 


(4.21) 


where  X.  are  the  eigenvalues  and  v.  the  corresponding  eigenvectors  of  M.  Thus, 

fji  I • 

if  v.  g is  not  small  for  small  eigenvalues  X. , the  step  in  the  direction  of  v.  is  large. 

1 tji  1 ^ 

Usually , for  small  eigenvalues , v.  g is  also  small , but  it  is  the  difference  of  large 
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and  almost  equal  numbers  . Therefore,  in  most  cases,  divergence  is  only  a numeri- 
cal problem.  Complications  arise  from  local  irregularities  in  the  likelihood  function. 
Several  remedies  can  be  used,  depending  on  the  problem  complexity: 


(a)  AO  = -aM  *g 

Let  A0'  = M_1g  (4.22) 

then  AO  = -ccAO*  (4.23) 

V 

Positive-definiteness  of  M could  be  lost  during  the  inversion.  This 
can  be  avoided  through  finding  A0  ' by  solving  the  following  linear 
equations 

MA0'  = g (4.24) 

The  accuracy  of  Eq.  (4.24)  also  depends  on  conditioning  of  M (ratio 
of  maximum  and  minimum  eigenvalues)  but  will  be  much  superior  as 
compared  to  Eq.  (4.22)  . 

(b)  Conditioning  of  M can  be  improved.  Premultiply  both  sides  of  Eq. 
(4.24)  by  A \ where 

A = diag  [v/M7r  • • • \/Mmm  ] (4-25) 

to  get 

(A_1MA-1)  (AAO  ')  = A_1g 
or 

M*A0*  = g*  (4.26) 
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All  elements  of  M*  are  between  -1  and  +1.  In  general,  this  matrix 
is  much  better  conditioned  than  M.  Equation  (4.26)  can  be  used 
to  solve  for  A0*  and  then  AG'  is  obtained  by  rescaling. 


(c)  The  rank  deficient  inverse  of  M can  also  be  used  . In  this  procedure , 
the  eigenvalues  of  the  matrix  M are  arranged  in  a decreasing  order  of 
magnitude  and  all  eigenvalues  below  a threshold  b are  neglected  in 
using  Eq.  (4.20)  for  the  inverse,  i.e.. 


A- 1 _ y 1 T 

M = 2j  -r—  V.V. 

. . A.  11 
1=1  1 


(4.27) 


where 


A.  > b 
i — 

i < k 

A.  < b 

i > k 

l 


(4.28) 


This  is  equivalent  to  searching  for  the  maximum  of  the  likelihood 
function  in  the  subspace  spanned  by  v.  (i  < k)  . If  this  process 
is  continued  until  the  end,  the  estimates  may  never  converge  to 
their  maximum  likelihood  values . It  should  also  be  noted  that  since 
only  large  eigenvalues  are  included,  step  A0  may  be  "small". 

(d)  Instead  of  neglecting  eigenvalues  of  M,  which  are  smaller  than  a 
certain  threshold  b,  these  smaller  eigenvalues  are  increased  to 
b.  Thus,  the  inverse  of  M is  approximated  as 


-1  1 T 1 m __ 

M = £ — v.v,  + - X v.vT 


i=l\  > * b 


i=k+l 


l i 


(4.29) 


In  this  case,  when  the  program  converges,  the  likelihood  function 
will  reach  a maxima.  Convergence  may  be  slow. 
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(e)  This  fix  of  the  Gauss-Newton  method  is  a two-step  procedure.  The 
first  step  is  taken  by  finding  M using  Eq.  (4.27).  Then,  another 
step  direction  is  computed  using  only  the  small  eigenvalues,  i.e. , 


A0 


2 


m 

-a  2 
i=k+l 


T 

v.v.  g 
l a b 


x: 


(4.30) 


Then  step  size  a is  chosen  by  a one  dimensional  search.  As  the 
algorithm  approaches  the  maximum,  it  may  be  worthwhile  to  carry 
out  a one  dimensional  search  in  each  v.(i  > k)  direction. 

The  last  approach,  when  used  with  the  normalized  information  matrix 
M*of  Eq.  (4.26),  is  one  of  the  best  and  is  recommended  if  the  number 
of  unknown  parameters  is  large.  In  most  cases  tried,  fix  (d)  was 
adequate . 

4.3.2  R elationship  Between  Subset-Regression  and  Rank-Deficient  Solutions 


The  subset-regression  method,  presented  in  Chapter  III  and  Appendix  B, 
and  the  rank  deficient  solution  of  Section  4.3,  are  closely  related  to  each  other. 

In  the  subset-regression  method,  an  equation  error  technique  is  used  in  conjunc- 
tion with  certain  statistical  tests  to  determine  the  set  of  parameters  which  signifi- 
cantly affect  the  input-output  relationship  for  the  given  data  and  the  class  of  models 
under  consideration.  In  other  words,  the  method  determines  which  parameters 
are  important.  Since  different  models  and  parameter  sets  may  be  adequate  for 
different  operating  regions,  the  importance  of  any  parameter  depends  upon  the 
input  and  output  data. 

The  rank  deficient  method  uses  input-output  data  and  the  model  structure 
to  determine  the  information  content  about  linear  combinations  of  parameters 
in  the  model.  It  works  in  parameter  space  finding  principal  directions  and 
information  about  the  linear  combination  of  parameters  representing  those  direc- 
tions (eigenvectors  and  eigenvalues  of  the  information  matrix).  When  the  rank 
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deficient  inverse  of  the  information  matrix  is  used,  the  parameter  subspace  with 
small  information  is  dropped  from  consideration.  The  optimization  procedure, 
then,  extremizes  the  likelihood  function  over  the  parameter  subspace  where 
enough  information  is  available  about  parameters  from  data.  This  could  be  used 
to  determine  the  model  structure  much  in  the  same  way  as  the  subset-regression 
technique.  However,  in  the  rank  deficient  method,  unlike  subset  regression, 
significant  linear  combinations  of  parameters  are  obtained.  In  the  numerical 
computation  procedure,  the  true  parameter  values  are  unknown  and,  therefore, 
only  the  conditional  information  matrix  for  assumed  parameter  values  is  available. 
This  information  matrix  could  be  considerably  different  from  the  true  information 
matrix  at  correct  parameter  values.  Thus,  dropping  a certain  parameter  subspace 
m earlV  stages  of  computation  may  lead  to  erroneous  conclusions.  This  factor 
reduces  the  applicability  of  this  method  in  model  specification,  since  the  optimal 
parameter  values  must  be  found  before  the  unimportant  parameter  directions 
can  be  determined . 


4.4  EXAMPLE  OF  THE  USE  OF  THE  MAXIMUM  LIKELIHOOD  IDENTIFICATION 

PROGRAM  (HIDENT) 

The  input  to  the  maximum  likelihood  program  are  noisy  measurements  and 
the  lequired  number  of  expansion  terms  in  all  aerodynamic  stability  and  control 
derivatives  . A priori  values  of  the  coefficients  in  each  of  these  expansion  terms, 
obtained  from  subset  regression  subroutine,  are  also  read.  This,  together  with 
the  equations  to  be  propagated  and  the  parameters  to  be  identified,  starts  the 
identification  subroutine,  denoted  as  INDENT. 

A typical  HIDENT  output  for  identifying  longitudinal  aerodynamic  coefficients 
is  presented  in  Tables  4.1  and  4.2.  The  stales  u,  w,  q,  and  0 are  propagated  by 
integrating  the  differential  equations  governing  them,  while  measurements  of  v, 
p,  r,  cp  and  vj/  are  assumed  to  be  the  values  of  the  corresponding  state  variables. 
The  measurements  of  a,  u,  z and  q are  used  to  identify  10  parameters  in  the 
longitudinal  equations  of  motion.  Many  other  parameters  have  nonzero  values 
but  are  not  identified.  For  example,  the  parameter  2 is  fixed  at  -5746  dur- 
ing identification  (i.e.,  is  not  identified).  a P 
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Table  4.2  is  a typical  iteration  step.  The  parameters  are  listed  and  then 
the  R matrix  (average  innovations  covariance)  is  given.  The  upper  half,  including 
the  diagonal,  is  the  R matrix  and  the  lower  half  gives  the  correlation  coefficients 
between  different  measurements.  Since  the  measurement  noise  covariance  is 
known  to  be  diagonal,  its  estimate  would  be 

R = -diagt  1 . 48  x l(f6,  3.97  x 10"5,  5.07  x l(f5,  8.58  x 10-6  ] (4.31) 

The  total  cost  for  these  sets  of  parameter  values  is  -15.5  and  represents  a 13%  cost 
decrease  from  the  parameters  in  the  last  step.  DJ  is  the  first  gradient  and  D2J  is 
the  second  gradient  (information  matrix)  of  the  cost  with  respect  to  parameters. 

The  eigenvalues  and  eigenvectors  of  the  information  matrix  are  given.  If  the 
range  of  eigenvalues  should  exceed  a ratio  of  10  ^ , one  of  the  remedies  suggested 
in  Section  4.3  is  used  depending  on  problem  complexity.  In  this  case,  since  the 

g 

range  of  eigenvalues  is  about  10  all  eigenvector  directions  are  included  in  the 
search.  The  inverse  of  the  information  matrix  (i.e.  , the  dispersion  matrix), 
the  new  parameters  and  the  parameter  step  are  shown  next.  The  terms  on  the 
diagonal  of  the  dispersion  matrix  give  parameter  error  variances,  the  square  root 
of  which  are  standard  deviations  of  estimates.  The  F-value  is  computed  as  the 
squared  ratio  of  the  parameter  estimate  to  its  standard  deviation. 

Thereafter,  we  proceed  to  the  next  iteration . This  process  continues  until 
convergence  occurs  (i.e. , the  change  in  cost  per  iteration  falls  below  a pre- 
specified tolerance  bound)  . 


4.5  SUMMARY 

The  maximum  likelihood  technique  of  parameter  identification  is  the  essential 
second  step  of  a method  for  quantifying  the  causes  of  aircraft  responses  in  the 
high  angle-of-attack  flight  regime.  Combining  the  two  advanced  techniques  of 
subset  regression  and  maximum  liklihood  is,  itself,  a complicated  procedure,  but 
one  which  is  necessary  for  isolation  of  the  phenomena  occurring  in  these  flight 
regimes.  As  discussed  in  Chapter  II,  however,  the  identification  algorithm  is  but 
one  part  of  the  entire  process.  For  this  effort,  the  consideration  of  input  design 
to  improve  parameter  estimates  has  also  been  important,  and  is  discussed  in  the 
following  chapter. 
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V.  INPUT  DESIGN 


The  principal  result  of  this  project  has  been  discussed  in  Chapter  III  (the 
model  determination  phase)  and  Chapter  IV  (the  maximum  likelihood  identification 
phase)  . Application  of  this  method  is  detailed  in  Chapter  VI.  Because  of  the  impor- 
tance of  control  input  design  in  achieving  those  results,  it  is  necessary  to  first 
review  various  aspects  of  input  specification  for  the  high  angle-of-attack  flight 
regime.  Requirements  for  such  input  designs  are  summarized  (Section  5.1)  and 
present  inputs  used  are  reviewed  (Section  5.2).  Three  approaches  for  high  angle- 
of-attack  input  design  are  discussed  (Section  5.3)  . These  are:  (1)  sequential 
linear  designs  for  locally  linear  approximations,  (2)  subset  regression  evaluation 
of  various  inputs,  and  (3)  analytical  nonlinear  input  designs.  Examples  of  pro- 
cedures (1)  and  (2)  are  given  (Section  5.4)  based  on  techniques  detailed  in 
Appendix  C.  The  third  technique,  developed  analytically  in  Appendix  C,  has  not 
been  applied  usefully  due  to  its  excessive  computational  requirements,  but  its 
basis  is  discussed  for  completeness . 

5 . 1 INPUT  REQUIREMENTS 


It  has  long  been  realized  that  the  ultimate  success  of  a flight  test  program 
depends  on  the  choice  of  inputs  used  to  excite  the  desired  motions  of  the  aircraft. 
Good  inputs  could  enhance  parameter  identifiability  and  improve  confidences  on 
parameter  estimates.  This  is  especially  true  for  parameter  identification  of  aircraft 
models  at  high  angles-of- attack . As  explained  in  Chapter  III,  the  aerodynamic 
derivatives  are  expanded  as  polynomials  in  angle-of-attack  and  sideslip  angle. 

This  leads  to  a large  number  of  unknown  parameters  . Inputs  must  be  chosen 
carefully  to  discriminate  among  different  parameters . There  are  many  other 
considerations  for  choosing  inputs  for  specific  flight  tests: 

Pilot  Acceptability:  The  pilot  should  be  able  to  safely  implement  these 
inputs  without  incurring  dangerous  aircraft  responses. 

Instrumentation:  The  inputs  should  consider  specific  instruments 
available,  their  dynamic  range  and  accuracy. 
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inputs  without  incurring  dangerous  aircraft  responses. 

b . Instrumentation:  The  inputs  should  consider  specific  instruments 
available,  their  dynamic  range  and  accuracy. 
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C.  Parameter  Identification  Procedure:  Certain  parameter  extraction 
methods  require  specific  inputs . There  is,  however,  no  constraint 
with  modern  powerful  procedures  like  the  maximum  liklihood  method 
used  with  the  modified  Gauss-Newton  optimization  technique. 

d.  Modeling  Assumptions:  The  input  should  be  tailored  to  the  model 

whose  parameters  are  to  be  identified.  In  high  angle-of-attack  flight 
regime,  the  aircraft  follows  different  models  for  different  regions  of 
operation.  The  input  should  ensure  that  the  system  operates  in  a 
region  where  the  model  is  valid . 


1 


e.  Objective  of  Parameter  Identification:  The  inputs  must  be  selected 

on  the  basis  of  the  ultimate  objective  of  parameter  identification;  e.g. , 
regulator  design,  response  prediction,  simulator  construction. 

f-  System  Integrity  Constraints:  The  aircraft  maneuvers  produced  by 
the  input  should  not  cross  the  design  stresses  of  various  aircraft 
components . 

g • Output  Sensitivity:  Measured  aircraft  response  resulting  from  the 

input  should  be  sensitive  to  the  parameters  of  interest  and  insensitive 
to  other  unknown  parameters  . 

5 . 2 CONVENTIONAL  INPUTS 


In  the  last  few  decades,  many  different  inputs  have  been  tried  for  aircraft 

parameter  identification.  Most  of  these  inputs  are  chosen  because  of  their  simplicity 

and  ability  to  yield  responses  from  which  the  relatively  small  number  of  basic 

linear  derivatives  (e.g.,  C , C , C0  , C , etc.)  may  be  estimated.  The 

m n 1 x n 

a q p r 

earliest  inputs  were  frequency  sweeps,  where  the  aircraft  is  excited  by  sinusoidal 


inputs  over  a range  of  frequencies  until  steady-state  is  reached  at  each  frequency. 
This  gives  the  frequency  response  from  which  parameters  of  a suitable  linear  model 
can  be  determined.  In  addition  to  being  very  time  consuming,  this  type  of  input 
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Recent  works  of  Gupta,  Mehra  and  Hall  [51-56]  have  given  considerable 
insight  into  the  nature  of  inputs  which  given  good  parameter  estimates  and  resolve 
identifiability  problems.  A discussion  of  techniques  to  determine  these  inputs  is 
given  in  Appendix  C.  Much  further  work  is  required  in  the  area  of  input  design 
for  linear  as  well  as  nonlinear  systems. 


5.3  INPUTS  FOR  HIGH  ANGLE- OF -ATTACK  FLIGHT  REGIME 


Stall 


5.3.1  Inputs  for  Identifying  Local  Parameter  Values 


ror  very  small  changes  in  a and  p during  a maneuver,  a model  assuming 
constant  values  of  aerodynamic  derivatives  is  usually  adequate,  as  mentioned 
in  Chapter  III.  The  techniques,  described  in  Appendix  C.2,  for  the  design  of 
optimal  inputs  for  linear  dynamic  systems,  can  be  used  to  specify  inputs  for 
identification  of  local  values  of  functions  representing  the  stability  and  control 
coefficients.  It  is  necessary  that  there  by  only  small  excursions  in  these  states 
which  have  maximum  effect  on  parameters  being  identified. 


5.3.2  Inputs  to  Identify  Nonlinear  Parameter  Variations 


1 


Many  me.thods  have  been  suggested  to  design  inputs  which  would  give 
good  identifiability  of  stability  and  control  coefficients  over  a range  of  states. 

One  way  could  be  to  carry  out  flight  tests  around  many  angles-of-attack  and 
sideslip  angles  covering  the  range  of  interest.  Inputs  based  on  linear  models 
are  designed  at  each  of  these  points,  as  explained  above,  lhis  would  give  an 
accurate  description  of  each  aerodynamic  derivative,  but  may  be  practically 
infeasible  because: 

1.  It  would  require  excessive  flight  testing  time  and  data  processing 
time. 

2.  It  would  require  carrying  out  flight  tests  in  regions  where  the  trim 
condition  cannot  be  reached  or  the  airplane  is  unstable  and/or  un- 
safe . 

A better  scheme  is  to  approximate  the  nonlinear  aerodynamic  derivatives 
by  polynomials  (or  some  other  truncated  series  of  complete  functions)  in  independ- 
ent variables.  Each  unknown  nonlinear  function  in  the  equations  of  motion  is 
thus  replaced  by  a set  of  parameters.  These  parameters  can  be  identified  in 
a variety  of  ways.  The  aerodynamic  derivatives  are  determined  from  a knowledge 
of  the  coefficients  of  terms  in  the  polynomial  expansions. 

For  small  ranges  of  a and  p , it  is  possible  to  determine  all  the  polynomial 
coefficients  in  one  single  experiment.  The  inputs  should  have  sufficient  amplitude 
so  that  tlie  motions  occur  in  more  or  less  the  entire  range  of  o and  P under  con- 
sideration. The  technique  for  determining  the  amplitude  is  given  in  Appendix 
C.4.2.1.  The  amplitude  depends  on  the  order  of  the  polynomial,  accuracy  desired, 
duration  of  the  experiment  and  a priori  estimate  of  the  polynomial  coefficients. 

To  determine  aerodynamic  derivatives  over  a wide  range  of  angle-of- attack 
(0<«<  30°)  and  sideslip  angle  (-20°  < p <20c),  it  is  necessary  to  perform  experi- 
ments starting  at  different  o,p  trim  conditions.  The  results  of  these  separate 
experiments  .ire  put  together  to  obtain  all  coefficients  in  the  polynomial.  The 
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algorithm  of  Appendix  C.4.2.2  finds  the  trim  conditions  and  the  duration  of  the 
experiment  at  each  trim  condition  to  produce  a good  estimate  of  the  nonlinear 
function  over  the  entire  range  of  interest.  It  is  possible  to  put  constraints  of  an 
infeasible  operation  over  a region.  The  technique  considers  regions  where  the 
stability  augmentation  system  (SAS)  is  required  for  pilot  safety  but  deteriorates 
parameter  estimates. 

A semi-empirical  technique  is  useful  in  the  design  of  input  signals  for 
identifying  parameters  in  complicated  systems  with  known  characteristics. 

As  a first  step,  analytic  input  designs  are  based  on  simplified  models,  using 
the  techniques  of  Appendix  C.  State  time  histories  are  generated  and  inputs 
are  evaluated  based  on  a more  accurate  but  complex  model.  This  would  give 
information  about  poorly  identifiable  directions  in  the  parameter  space  and  about 
poorly  excited  modes.  This  is  usually  done  using  simple  programs  which  are 
not  necessarily  recommended  for  data  reduction  from  actual  flight  tests . It  is  to  be 
noted  that  input  evaluation  can  be  done  without  generating  simulation  data  and 
going  through  the  identification  process.  In  the  implementation  used  here,  the 
optimal  subset  regression  program  is  used  to  determine  identifiable  directions 
resulting  from  a certain  input.  A knowledge  of  the  deficiencies  in  the  chosen 
input,  together  with  known  system  behavior,  is  used  to  modify  the  input.  The 
process  is  repeated  until  an  acceptable  input  is  obtained. 

Some  analytical  results  for  the  design  of  inputs  for  nonlinear  systems  are 
presented  in  Appendix  C.  Many  of  them  lead  to  complicated  problems  which  are 
difficult,  if  not  impossible,  to  solve  with  present  computational  facilities.  They 
do  give  an  excellent  insight  into  the  nature  of  inputs  which  could  give  good 
identifi ability  of  unknown  parameters.  If  the  form  of  a nonlinear  stability  or 
control  derivative  is  now  known  a priori,  it  is  necessary  to  operate  the  system 
over  the  entire  range  of  angle-of-attack  and  sideslip  angle  of  interest.  In  other 
words , extrapolation  could  give  misleading  results  when  a nonlinear  function  is 
approximated  by  a finite  series  of  complete  functions  over  a range  of  the  independ- 
ent variables.  Secondly,  most  input  energy  should  be  concentrated  over  the  test 

rc,  ime . If  the  local  variation  of  C with  angle-of-attack  is  to  be  studied,  it  is 
° m 

not  good  to  use  inputs  which  give  large  excursions  in  angle-of-attack,  Thl*.  will 
' c demonstrated  through  examples  in  the  next  section. 
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5.4  EXAMPLES  OF  INPUT  DESIGN  IN  HIGH  ANGLE -OF- ATTACK  FLIGHT 
REGIME 


5.4.1  Static  Pitching  Moment  Characteristics 


One  of  the  most  important  phenomena  of  the  stall/post-stall  regime  is 
pitch-up".  Its  principal  cause  is  the  positive  C which  occurs  close  to  stall  [1] 


Unlike  most  lateral  aerodynamic  nonlinearities,  the  affected  mode  is  the  rela- 
tively simple  short  period  response  which  can  be  conceptually  treated  as  a second 
order  system  whose  stiffness  is  controlled  by  C 


The  primary  objective  is  to  investigate  the  following  questions 


What  flight  test  conditions  are  required  to  determine  the  presence 
of  C >0? 


What  order  polynomial  is  required  to  represent  the  versus  a 
characteristic  over  the  test  regime? 


In  order  to  resolve  these  issues,  several  initial  conditions  and  inputs  were 
applied  to  the  F-4  simulation;  the  subsequent  responses  were  passed  through  the 
subset  regression  program.  Since  the  principal  effect,  of  the  aerodynamic  non- 
linearity was  known  to  be  on  pitch  (short  period)  response,  only  ELQ,  the  pitch 
acceleration  regression  variable,  was  used. 


Initial  conditions  and  inputs  used  were  the  following 


Trim  at  1 ,10  ,20  angle-of~ attack , doublet  stabilator  1 
5°,  respectively. 


Trim  at  10  angle- of- attack 


Trim  at  13  angle-of-attack , sinusoidal  stabilator  (5 
amplitudes)  . 


The  results  of  perturbation  inputs  in  anticipated  linear  regimes  were  developed 
first  [1]  . These  regimes  were  low  angle-of- attack  (a  = 1°)  , pre-stall  (a  = 10°) , 
and  post-stall  (a  = 20°)  . The  data  from  these  maneuvers  was  passed  to  the  sub- 
set regression  program.  The  resulting  coefficient  model  was  evaluated  on  the 
basis  of  the  quantitative  "fit"  of  the  response  from  the  estimated  model  to  the 
actual  simulation  response.  The  results  are  shown  in  Table  5.1. 


The  next  series  of  longitudinal  inputs  [2]  were  made  with  a stabilator  ramp 
of  l°/sec  to  slowly  push  the  aircraft  from  the  linear  regime  through  the  known 
nonlinear-  Cm  versus  a characteristics  at  about  15°  (starting  from  10°).  In 
addition,  sinusoidal  inputs  were  used  (case  (3))  . 


The  basic  result  of  these  runs  was  that  the  C versus  a nonlinear  character 

m 

istics  required  data  replication  for  successful  determination.  Only  the  last  series 
of  runs,  (3),  provided  this  replication.  The  physical  reason  for  this  requirement 
is  based  on  the  negative  spring  interpretation  of  the  static  moment  characteristic. 


TABLE  5 .1 

Stabilator  Doublet  (These  results  are  a compilation  of  several  experiments) 


ANGLE  OF 
ATTACK 


MODEL 


+1  stabilator  doublet,  7 second  data  length 


C +C  a+C  6 +C  q matched  (virtually)  100%  oi 
m m mc  s mn  ..  . . / ... 

o a o q pitching  moment  variation 


+ 1 stabilator  doublet 


8 +C  q matched  99 . 8%  of  pitching 
q moment  variation 


+5  stabilator  doublet,  7 second  data  length 


matched  99. 1%  of  pitching 
moment  variation 


If  the  initial  rates  and  accelerations  with  which  the  aircraft  enters  the  unstable 
region  are  too  low,  the  response  to  a doublet  or  step  would  be  unstable,  forcing 
the  response  to  a stable  region.  If  the  initial  accelerations  are  too  high,  the 
aircraft  inertia  overwhelms  the  static  moment  and  the  data  doesn't  reflect  the 
nonlinearity.  Although  several  sets  of  initial  conditions  and  time-varying  inputs 
may  exist  for  isolating  the  characteristics,  no  extensive  experimentation  was 
performed  with  other  than  a sinusoidal  stabilator  input. 

The  results  for  the  sine  wave  inputs  are  shown  in  Figures  5.1  through  5.4. 
Figure  5.1  shows  the  10°  amplitude  sine  wave  input  and  the  resulting  aircraft 
response  (a,  q,  ELQ)  . This  input  was  based  on  a trim  angle- of- attack  of  13.5°. 

It  is  recognized  that  such  a trim  may  be  difficult  to  establish  in  a flight  test  since 
the  aircraft  is  close  to 

max 

The  data  from  this  run  was  passed  through  the  regression  problem  for  two 

maximum  allowable  polynomials — a fourth  order  and  a ninth  order.  The  results 

are  shown  in  Figure  5.2.  Compared  with  the  (known)  simulation  C versus  a, 

m 

this  input  did  not  come  acceptably  close  to  the  "true"  value. 

Based  on  the  mechanism  discussed  above,  a new  input  with  half  the  ampli- 
tude and  frequency  was  applied,  the  response  for  which  is  shown  in  Figure  5.3. 
This  input  should  provide  a more  acceptable  regression  result,  which  is  evident 
in  the  new  subharmonic  response  in  ELQ.  Passing  the  data  of  this  new  input 
through  the  regression  program  (again  with  two  polynomial  possibilities)  produced 
the  result  shown  in  Figure  5.4.  Clearly,  a higher  order  polynomial  can  be  identi- 
fied (locally)  using  this  input. 

5.4.2  Yaw  Stiffness  Variation  with  Angle-of- Attack 

Yawing  moment  produced  by  sideslip  angle  (C  ) is  an  important  parameter 

"P 

in  determining  stability  of  the  lateral  modes  of  the  aircraft.  A positive  C usually 

np 

implies  a stable  dutch  roll  mode  and  vice  versa.  C varies  with  angle-of- attack . 

np 

Figure  5.5  shows  a plot  of  C as  a function  of  angle-of-attack  a at  zero  sideslip 

np 

angle.  It  is  generally  positive  in  the  low  angle-of-attack  flight  regime  and  changes 
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sign  beyond  the  stall  angle.  Also  shown  in  Figure  5.5  is  a least  squares  fourth 
order  polynomial  fit  to  the  simulation  data. 


I 


A computer  program*  has  been  written  to  determine  the  angles-of- attack  at 
which  C must  be  identified  locally  to  give  a good  estimate  of  the  C curve  for 

o "P  o ^ 

0 < a < 30  . It  is  assumed  that  operation  is  possible  throughout  this  range  of 

angle-of- attack.  If  the  angle-of-attack  is  greater  than  20° , however,  it  is  neces- 
sary to  use  the  stability  augmentation  system  to  make  the  airplane  stable  enough. 
This  use  of  SAS  at  high  angles-of- attack  results  in  poorer  parameter  estimates 
for  the  same  duration  of  the  experiment  time.  The  information  available  about 
C at  any  a is  assumed  to  increase  linearly  with  the  experiment  time  around  that  a. 

P 

The  use  of  SAS  increases  the  standard  deviations  by  a factor  of  three  for  the  same 

experiment  time.  The  information  about  C for  a unit  experiment  time  at  different 

np 

angles-of-attack  is  shown  in  Figure  5.6. 

To  identify  coefficients  in  the  polynomial  which  approximate  C as  accu- 

“P 

rately  as  possible,  a 100  sec.  long  experiment  time  should  be  divided  as  shown 
in  Table  5.2.  Many  different  designs  would  give  the  same  information  matrix 
for  unknown  parameters.  The  design  with  a minimum  number  of  points  where 
the  experiment  need  be  carried  out  is  chosen.  Notice  that  to  identify  coefficients 
in  an  nth-order  polynomial,  it  is  necessary  to  find  the  value  of  the  polynomial 
at  no  less  than  (n+1)  points.  Therefore,  five  is  the  minimum  number  of  points 
where  the  experiment  need  be  carried  out. 


5.5  SUMMARY 

The  role  of  input  design  for  parameter  extraction  is  significant  for  high  angle- 
of-attack  applications.  Linear  input  designs  are  useful,  particularly  for  local 
linear  flight  regimes  and  for  determining  flight  test  schedules  that  minimize  time 
to  obtain  estimates  with  such  linear  inputs.  For  highly  nonlinear  regimes,  however, 
iterative  use  of  the  model  determination  program  will  aid  in  specifying  required 
inputs.  Analytical  design  methods  for  nonlinear  regimes  can  be  formulated,  but 
appear  to  present  computational  problems  which  must  be  resolved  by  further  study. 

“ 

This  is  a modification  of  a program  written  by  SCI  for  NASA  Edwards  Flight 
Research  Center,  Edwards,  California  under  Contract  NAS  4-2068. 
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VI.  IDENTIFICATION  IN  THE  STALL/POST-STALL 
HIGH  ANGLE-OF-ATTACK  REGIME 


6.1  INTRODUCTION 

In  this  chapter,  evaluation  of  the  integrated  j)arameter  identification  process 
is  carried  out  using  simulated  aircraft  flight  test  data.  The  simulation  employed 
is  that  of  the  F-4  aircraft  described  in  ref.  [1]  , Chapter  II  and  Appendix  D. 

State  and  measurement  model  equations  used  in  this  maximum  likelihood  procedure 
are  detailed  in  Appendix  E.  Nine  state  equations  and  eight  measurement  equations 
containing  thirty-three  nonlinear,  aerodynamic  coefficients  compose  the  identifica- 
tion model.  (Each  of  these  coefficients  are  further  expanded  as  discussed  in  Chapters 
III  and  IV  . ) 

In  the  identification  procedure,  the  subset  regression  program  is  used 
on  the  simulated  flight  test  data  to  identify  the  model  structure  and  give  initial 
estimates  of  the  selected  aerodynamic  coefficient  polynomial  expansion  parameters. 
The  maximum  likelihood  program  is  then  used  to  refine  the  parameter  estimates 
obtained  from  the  regression  analysis  to  yield  the  final  parameter  estimates  and 
a measure  of  the  confidence  associated  with  those  estimates.  The  detailed  flow 
chart  of  the  entire  process  is  illustrated  in  Figure  6.1. 

In  addition  to  evaluating  the  identification  pi*ocedure  on  a variety  of  simulated 
flight  test  experiments,  investigations  are  also  carried  out  to  assess  the  effects 
of  the  following: 


Different  levels  of  measurement  noise 

Process  noise 

Control  input  variations 

Identified  parameter  set  size 

Data  length 


The  first  two  items  illustrate  the  robustness  of  the  procedure  to  corrupting  influ- 
ences. In  general,  increased  noise  levels  degrade  both  the  final  estimates  of 
parameters  and  the  ability  of  the  regression  to  identify  model  structure.  The 
last  three  items  show  the  intrinsic  limitations  imposed  on  parameter  estimates 
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INTEGRATED  PARAMETER  IDENTIFICATION  PROCEDURE 
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Figure  6. 1 Operational  Flow  Chart  of  Integrated  Parameter 
Identification  Process 


by  the  information  content  of  the  data.  Increasing  the  information  content  by 
appropriate  input  design  and/or  increasing  the  amount  of  data  in  time  or  numbers 
of  measurements  will  increase  the  number  of  parameters  that  can  be  identified 
with  a significant  level  of  confidence. 

6.2  EVALUATION  OF  THE  INTEGRATED  IDENTIFICATION  PROCESS 

The  objective  of  the  integrated  parameter  identification  procedure  is  to 
determine  values  for  unknown  or  poorly  known  parameters  of  a mathematical 
model  of  a process  by  matching  data  from  the  real  process  with  the  model's 
predicted  results.  For  the  process  considered  in  this  work,  the  parameters 
are  the  polynomial  expansion  terms  in  a and  p of  the  aerodynamic  coefficients 
in  the  aircraft  rigid  body  dynamics  model.  Aircraft  dynamic  data  is  obtained 
by  means  of  a simulation  using  a much  more  detailed  model. 

6.2.1  Selection  of  Test  Conditions 

The  identifi ability  of  the  stall/post-stall  regime  is  the  prime  objective  of 
the  applications  conducted  for  the  identification  process.  This  angle-of-attack 
regime  (10°  to  25°  and  beyond)  is  characterized  by  multiple  nonlinearities  in  the 
pitch,  roll,  and  yaw  moments  and  forces.  The  central  objective  of  this  applica- 
tion is  that  resolution  of  the  identifi  ability  problems  for  this  range  is  the  prime 
requisite  for  any  high  angle-of-attack  identification  procedure. 

The  evaluation  procedure  which  has  been  determined  is  based  on  the 
following  considerations: 

1.  Selection  of  Inputs:  Not  all  inputs  will  sufficiently  excite  the  air- 

craft to  induce  nonlinear  forces  and  moments . These  results  were 
noted  in  Chapter  III.  A standard  set  of  inputs  is  selected  for  this 
purpose . 

2.  Selection  of  Data  Length:  In  general,  the  longer  the  data  length  , 

the  better  will  be  the  identification  accuracy  of  those  parameters 
which  can  be  identified.  At  high  angl  es-of-  attack , however,  the 
amplitude  of  responses  may  prohibit  extensive  time  at  a particul 


! 
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flight  condition.  This  consideration  led  to  a data  length  of  10  seconds, 
with  a sample  rate  of  10/second. 


Selection  of  Primary  Coefficients:  The  coefficients  which  affect 

aircraft  response  may  be  classified  as  primary  or  secondary.  Primary 
coefficients  are  those  which  most  affect  aircraft  response  and  which 
are  of  greatest  interest  to  the  test  engineers.  Secondary  coefficients 
are  those  which  have,  in  general,  smaller  effects  on  aircraft  response, 
but  which,  for  certain  requirements  (e.g.  , SAS  or  handling  quality 
studies),  may  be  of  interest.  For  the  applications  discussed  here, 
the  static  force  and  moment  coefficients  are  considered  primary 


(e.g.,  C , C , C , C , C , Cp).  The  subset  regression  program 


x'  y'  z'  m'  n' 
is  the  essential  step  in  the  definition  of  primary  or  secondary  coefficients. 


Tabular  summaries  of  the  selected  test  conditions,  inputs  and  noise  levels 
are  given  in  Figure  6.2.  Three  basic  classes  of  test  condition  are  given — lateral 
mode  excitation,  longitudinal  mode  excitation,  and  coupled  longitudinal-lateral 
mode  excitation. 


In  the  first  class,  control  inputs  are  chosen  to  excite  the  lateral  direc- 
tional modes  of  the  aircraft: 


• Lateral  velocity , v 

• Roll  rate,  p 

• Yaw  rate,  r 


Applied  controls  are  a half  period  sine  wave  pulse  in  rudder,  8^,  of  10°  amplitude 
and  period  duration  of  2.5  seconds  followed  by  a full  period  sine  wave  doublet 
in  ailerons  of  7°  amplitude  and  period  duration  of  2. 5 seconds.  Total  data  length 
is  10  seconds  real  time,  sampled  every  0. 1 seconds  to  give  100  points  for  each 
measured  variable.  Initial  flight  conditions  are  specified  in  Figure  6.2  under 
Case  1.  Cases  2 and  3 use  the  same  control  sequence  as  Case  1,  but  with  different 
noise  characteristics.  Case  4 adds  a simultaneously  applied  stabilator  input 
to  force  the  aircraft  through  a larger  range  of  angle-of- attack . 
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INITIAL  CONDITIONS  CONTROL  INPUTS  (from  initial  conditions)  I NOISE 


Lateral  Response  Test  Cases 


INITIAL  CONDITIONS  CONTROL  INPUTS  (from  initial  conditions)  NOISE 


Figure  6.2c  Longitudinal/Lateral  Test  Cases 


— 


The  second  class  primarily  involves  the  longitudinal  mode  with  an  input 
chosen  to  excite: 

• Pitch  rate , q 

• Vertical  velocity,  w 

• Longitudinal  velocity , u 

Stabilator  (elevator)  input  (Figure  6.2b)  is  a continuous  sine  wave  superimposed 
on  the  initial  constant  stabilator  angle.  With  a period  of  4 seconds,  the  input 
is  designed  to  be  approximately  at  the  short  period  mode  frequency  of  the  aircraft 
and,  therefore,  yield  the  most  possible  information  about  the  parameters  governing 
that  mode.  Amplitude  of  the  input  is  -4°.  Again,  the  data  length  is  10  seconds 
real  time  and  there  are  100  sample  points.  Table  6.2  contains  a complete  description 
of  this  experiment  under  longitudinal  case  1.  Case  2 is  also  a sinusoidal  stabilator 
input,  but  with  smaller  amplitude. 

Finally,  inputs  to  excite  more  coupled  responses  are  specified  as  longitudinal- 
lateral  cases  1 and  2.  These  inputs  are  shewn  in  Figure  6.2c. 

Lateral  case  1 and  longitudinal  case  1 are  "baseline"  cases  which  are  evalu- 
ated first  to  establish  a reference  for  comparison  with  input,  noise  and  data  length 
effects . 

These  tests  are  performed  in  the  stall/post-stall  regime  so  that  nonlinear 
a and  p expansion  terms  in  some  of  the  aerodynamic  coefficients  should  be  required 
for  large  enough  excursions  in  a and  p.  To  assure  coefficient  a and  p variation, 
relatively  large  control  inputs  are  necessary.  An  important  point  is  that  it  is 
difficult  to  accurately  identify  parameters  affecting  specific  aircraft  motion  and 
characteristics  if  the  data  do  not  contain  enough  information  about  those  motions 
and  characteristics.  Hence,  the  nonlinear  regimes  must  be  entered  (which 
usually  implies  large  control  inputs)  in  order  to  identify  nonlinear  effects. 

6.2.2  Application  of  Subset  Regression 

Univariate  subset  regression  is  used  with  the  simulated  data  in  the  para- 
meter identification  procedure  to  develop  a reduced  set  of  significant  parameters 
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from  the  full  set  allowed  for  identification.  Initial  estimates  of  the  significant 
parameters,  although  biased  due  to  measurement  noise,  are  also  determined. 

The  regression  formulates  an  identification  model  by  assigning  estimated  values 
of  zero  to  parameters  which  have  little  or  no  effect  in  matching  the  given  data. 

As  discussed  in  Chapter  III,  subset  regression  is,  therefore,  essentially  equivalent 
to  a rank  deficient  least  squares  equation  error  identification  procedure  applied 
to  a single  equation. 

6 . 2 . 2 . 1 Independent-Dependent  Variable  Determination 

For  aircraft  identification,  the  regression  program  is  applied  successively 
to  each  of  the  six  differential  equations  described  in  Chapter  III  modeling  motion 
in  each  of  the  six  degrees  of  freedom.  Those  equations  are  of  the  form: 


u = fj(x,0) 

(a) 

v = f2(x,0) 

(b) 

% 

W = f3(x,0) 

(c) 

P = 

(d) 

q = f5(x,0  ) 

(e) 

r = f6(.x.,0) 

(f) 

where 


x = col{u,v,w,p,q,r,<p,0,\|r}  * (6.2) 

and  the  Euler  angles  given  in  Appendix  E. 

In  order  to  calculate  the  unknown  parameters,  0,  using  regression,  Eqs. 

(6. 1)  and  (6.2)  must  be  put  into  the  form 

* 

Note  that  0 is  the  vector  of  unknown  parameters,  while  the  scalar  0 is  the  aircraft 
pitch  angle . 


" ' " "■  ' ' ; ' ■ 


Yi  = ® i + ec 


(6.3) 


where  the  independent  variables  are  modified  versions  of  x^,  and  £ is  a vector 
of  linear  and  nonlinear  dependent  variables  composed  from  the  state  variable, 
x,  subset  u,  v,  w,  p,  q,  and  r.  8^  is  an  unknown  but  constant  parameter. 

The  transformations  to  produce  the  variables  in  Eq.  (6.4)  are  given  in  Chapter  III. 


The  measurements  available  from  the  simulated  flight  data  are  (ref. 
Appendix  D): 
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These  measurements  coupled  with  the  assumed  known  values  for 

• aircraft  mass  and  moments  of  inertia, 

• engine  thrust  and  moments, 

• aerodynamic  constants,  air  density  and  velocity, 

are  sufficient  to  calculate  measured  values  for  all  independent  and  dependent 
variables  in  Eq.  (6.2)  according  to  the  equations  of  Appendix  E.  There  are 
six  independent  variables  y corresponding  to  six  left  hand  side  variables  of 
the  model  (c.f.  eqs.  of  Chapter  III)  . The  dependent  variable  list  has  30  elements. 
As  Table  6.1  indicates,  however,  not  all  30  variables  are  admissable  for  regression 
against  each  independent  variable. 

The  dependent  variables  involving  a and  P are  transformed  to  variables 
(u-a)  and  (P“P),  where  a and  (T  are  average  values  for  the  particular  experiment 
being  evaluated.  This  procedure  "centers"  the  data  for  a and  p variations  and 
is  helpful  in  reducing  correlations  between  dependent  variables  that  are  functions 
of  a and  p.  High  correlations  between  dependent  variables  make  it  difficult 
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for  the  regression  to  distinguish  significant  differences  between  dependent  variables. 
The  result  is,  for  variables  that  are  powers  of  a and/or  p,  that  high  power  variables 
may  be  substituted  for  low  power  variables  for  statistically  small  differences 
in  significance  levels.  Such  a situation  is  undesirable  since  the  low  power  variable 
is  usually  preferable  for  reasons  such  as  the  increased  error  of  high  power  terms 
due  to  noise. 

The  stall/post-stall  regime,  however,  may  require  such  high  order  terms 
simply  to  describe  the  nonlinearities.  It  is  the  tradeoff  between  the  need  to  account 
for  high  order  nonlinear  contributions , and  the  desire  to  use  the  lowest  order 
terms  possible,  which  constitutes  a major  problem  with  any  model  estimation 
technique . 

6. 2. 2. 2 Evaluation  of  Regression  Results 

Univariate  subset  regression  provides  several  statistical  controls  for  deter- 
mination of  the  most  significant  variables  for  inclusion  in  the  fit  of  a particular 
dependent  variable: 

a.  Error  test  on  overall  fit 

b.  F-value  for  inclusion  of  a variable 

c.  F-value  for  deletion  of  a variable 

d.  Tolerance  or  correlation  of  variables 

The  error  test  is  a measure  of  fit  that  varies  from  0 for  worst  to  1.0  for  best  fit. 
F-value  is  a statistic  calculated  for  each  variable's  correlation  coefficients  to  deter- 
mine  the  relative  significance  and  confidence  of  a coefficient  value.  Specifically, 
the  F-value  is  the  ratio  of  a numerical  parameter  mean  divided  by  the  standard 
deviation  of  that  estimate.  A high  F-value  implies  a low  coefficient  standard  devia- 
tion and,  hence,  good  confidence  in  the  coefficient  estimate.  Tables  of  F-distribu- 
tions  are  used  to  choose  a critical  F-value  for  a particular  confidence  level  to 
include  or  delete  a variable  in  the  fit  given  the  number  of  data  points  to  be  regressed 
and  the  number  of  admissible  variables.  Finally,  tolerance  is  a measure  of  the 
maximal  correlation  of  variables  not  included  in  the  fit  variables  or  combinations 
of  variables  already  included  in  the  fit.  Tolerance  levels  vary  between  1 .0  for 
perfect  correlation  to  0 for  no  correlation.  A more  detailed  explanation  of  these 
statistics  is  given  in  Appendix  B . 

\ 
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For  the  test  conditions  described  in  Section  6.2,  the  various  regression 
control  parameters  are  given  in  Table  6.2.  The  regression  procedure  adds  variables 
to  the  fit  one  at  a time  in  decreasing  order  of  significance  until  either  the  F value 
or  tolerance  of  unincluded  variables  is  below  the  preset  values , or  until  the  fit 
error  test  is  above  another  preset  value.  The  purpose  is  to  identify  as  many 
significant  aerodynamic  parameters  as  possible  and  not  to  merely  match  a time 
history.  Hence,  this  error  test  criteria  is  set  very  high  in  order  that  all  variables 
will  be  included  which  meet  the  99%  coefficient  confidence  level  implicit  in  the 
F value  limit  and  the  tolerance  criteria. 

TABLE  6 .2 

Regression  Statistical  Control  Values 


Error  Test  Criteria 

0.999 

F-value  for  inclusion  of  a variable 

3.0 

(99%  confidence  level) 

F-value  for  deletion  of  a variable 

3.0 

(99%  confidence  level) 

Tolerance 

0.01 

Number  of  data  points 

100 

Number  of  admissible)  Lateral 

26 

dependent  variables)  Longitudinal 

22 

The  desire  to  determine  as  many  significant  variable  coefficients  as  possible 
must  be  tempered  to  some  extent  in  that  addition  of  the  last  few  variables  may 
only  contribute  marginally  to  the  fit.  Statistically  these  marginal  fit  improve- 
ments may  not  be  significant,  and,  hence,  the  variable  should  not  be  included. 

The  justification  for  this  is  that  although  the  fit  to  the  available  data  is  improved, 
the  number  of  degrees-of-freedom  in  the  fit  relative  to  the  degrees-of-freedom 
in  the  data  is  increased  and,  therefore,  the  confidence  bands  for  predictions 
to  match  new  data  could  be  reduced.  The  F-ratio  for  the  overall  fit  (see  Appendix 
B)  measures  this  tradeoff  between  fit  quality  and  degrees-of-freedom  in  the  fit 
versus  degrees-of-freedom  in  the  data.  Therefore,  in  selecting  the  best  variable 
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set  from  the  regression  for  use  in  the  maximum  likelihood  procedure,  the  maximal 
F-ratio  value  determines  which  of  the  sequential  fit  attempts  produced  by  the 
regression  should  ultimately  be  used.  Figures  3.4  and  3.5  of  Chapter  III  show 
the  cutoff  criterion  for  selecting  parameters  to  be  identified.  In  general,  all 
parameters  added  to  the  regression  up  to  and  including  the  maximum  F-value 
are  identified  by  the  maximum  likelihood  algorithm  (c.f.  Figs.  3.4  and  3.5). 

In  summation,  the  criteria  for  using  the  stepwise  regression  are: 

1.  Set  F-level  limit  for  variable  inclusion/deletion  to  the  99%  confidence 
level  based  on  available  data  and  admissible  variable  set. 

2.  Set  tolerance  limit  at  about  0.01. 

3.  Set  error  test  criteria  such  that  the  regression  steps  terminate  only 
on  F-level  or  tolerance  criteria. 

4.  From  among  the  regression  steps  performed,  use  results  from  that 
step  which  has  the  maximal  F-ratio  for  the  overall  fit. 

The  coefficient  values  for  the  regression  variables  obtained  according  to  the 
above  criteria  are  then  used  as  a priori  parameter  estimates  in  the  maximum  likeli- 
hood procedure. 

Regression  results  for  lateral  case  1 and  longitudinal  case  1 are  summarized 
in  Table  6.3.  In  the  longitudinal  case,  note  the  relatively  few  variables  and 
their  low  F-values  in  the  u equation  fit.  The  low  confidence  is  a result  of  the 
low  information  present  in  the  data  interval  which  is  short  relative  to  the  phugoid 
mode  time  constants.  By  contrast,  in  the  lateral  example,  the  v equation  fit  also 
has  only  a few  variables,  but  the  high  F-value  indicates  good  confidence  in  the 
coefficient  values.  The  conclusion  is  that  the  motion  in  the  data  is  well  described 
by  just  the  few  linear  terms  selected  by  the  regression  procedure. 
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TABLE  6.3 

Regression  Results  for  Demonstration  Cases 
Lateral  Response  Case  1 


I'AltAMPTEU 


I'AllAMKTKR 

KNTJMATK 

‘5.06  x 

ID' 1 

3.69  x 

10-4 

-6.97  x 

10"2 

-3.01  x 

rvj 

i 

o 

2.45  x 

101 

1.05  x 

10"1 

-1.69  x 

101 

4.95  x 

10'1 

1.51 

-2.47  x 

10"1 

2.79  x 

10"1 

1.41  x 

101 

3.56  x 

io'3 

-3.28 

-2.52  x 

10"2 

2.02  x 

io'3 

-4.18  x 

10‘2 

4.64  x 

IO-2 

1.26  x 

io-2 

2.47 
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4.19  x 
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-2.66  x 

101 

-3.31  x 
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TABLE  6.3  (Continued) 
Regression  Results  for  Demonstration  Cases 

Longitudinal  Response  Case  1 


EQUATION 


OVERALL  FIT 
F -RATIO 


PARAMETER 


PARAMETER 

ESTIMATE 


8.56  x 10 


-1.49  x 10 


-9.56  x 101 


2.01  x 10J 


-2.13  x 10 


-7.60  x 10 


-4.53  v 1(1 


-8.40  x 10‘ 


-8.49  x 10' 


-5.06  x 10 


-1.48  x 10‘ 


-4.42  x 10 


PARAMETER 

F-VALUE 


5.56  x 102 
8.60  x 101 


F > 4%  F 


The  other  equation  fits  contain  more  variables,  including  several  nonlinear 
polynomial  expansion  terms.  High  confidence  coefficient  estimates,  however, 
are  primarily  associated  with  linear  variables,  indicating  that  even  for  the  fairly 
large  ct,p  excursions  (see  Figure  6.2)  encountered  in  these  examples,  linear 


9 


h 


. 

effects  predominate.  It  is  expected  that  if  the  operating  range  is  reduced  enough, 
the  aircraft  motion  should  be  described  quite  well  by  the  linearized  versions 
of  the  aerodynamic  coefficients  in  the  model  equations . Nonlinear  effects  should 
basically  be  second-  and  lower-order  type  effects  superimposed  on  the  basic 
model  linearized  about  the  operating  point. 

6.2.3  Maximum  Likelihood  Application 

Parameter  estimation  by  the  maximum  likelihood  method  in  the  output  error 
formulation  (see  Chapter  IV)  is  the  final  step  in  the  identification  procedure 
used.  Maximum  likelihood  estimation  provides  refined,  unbiased  estimates  of 
the  parameter  set  defined  by  the  regression  results.  Actually,  the  maximum  likeli- 
hood portion  of  the  procedure  could  be  applied  without  the  regression,  provided 
reasonable  a priori  estimates  of  the  parameters  could  be  obtained  f *om  another 
source,  such  as  from  wind  tunnel  data.  Even  if  such  a priori  estimates  are  available, 
however,  there  is  a penalty  in  terms  of  required  computation.  Because  a priori 
knowledge  of  the  subset  of  significant  parameters  for  the  data  being  evaluated 
is  not  available,  estimates  must  be  attempted  for  the  full  parameter  set.  Since 
computation  to  obtain  parameter  estimates  by  the  maximum  likelihood  function 
optimization  is  more  lengthy  than  for  the  regression,  it  is  more  efficient  to  develop 
initial,  albeit  biased  and  less  confident,  estimates  and  significance  ranking  by 
the  regressions  and  then  use  the  results  as  a starting  point  for  the  maximum  likeli- 
hood procedure. 

6.2.3. 1 Parameter  Set  Decomposition  for  Maximum  Likelihood  Identification 

As  discussed  above,  the  computational  requirements  of  maximum  likelihood 
make  it  very  desirable  to  concentrate  as  much  as  possible  on  estimates  of  only 
those  parameters  which  can  be  estimated  well  from  the  available  data.  Since 
the  computation  requirements  are  roughly  proportional  to  both  the  number  of 
parameters  and  the  number  of  model  states  propagated,  estimating  parameters 
which  are  ill-specified  by  the  data  materially  increases  computation  with  only, 
at  best,  marginal  increases  in  useful  results.  In  fact,  as  Section  6.4  will  illustrate, 
over-parameterization  can  actually  decrease  the  quality  of  all  the  estimates. 
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First  of  all,  many  parameters  can  be  isolated  by  examination  of  the  data. 

If  it  is  apparent  that  there  is  only  minor  excitation  and,  hence,  little  information 
corresponding  to  a particular  identification  model  state,  the  state  equation  for 
that  model  and  its  associated  parameters  can  be  removed  from  the  identification 
procedure.  Either  a constant  or  a measured  time  history  of  that  state  can  then 
be  substituted  in  the  remaining  state  equations  where  there  is  cross-coupling 
to  the  removed  state.  This  technique  can  be  applied  in  much  parameter  esti- 
mation work  involving  aircraft  where  the  longitudinal  and  lateral  degrees  of 
freedom  can  frequently  be  decoupled  for  inputs  confined  to  one  plane  or  the  other. 
The  two  example  experiments,  lateral  case  1 and  lateral  case  2,  illustrate  this 
state  equation  elimination,  where  in  place  of  the  original  nine  state  equations 
for  the  complete  identification  model,  only  four  equations  are  required  in  the 
longitudinal  experiment  and  five  in  the  lateral.  The  cross- coupling  to  eliminated 
states  is  satisfied  by  using  measured  time  histories  of  the  states  nc.t  propogated . 

Isolation  of  a reduced  set  of  significant  parameters  by  equation  elimination 
can  generally  be  applied  before  using  subset  regression  on  the  data.  After 
regressing  the  data,  if  the  results  for  a particular  equation  reveal  a low  F-ratio 
for  the  overall  fit  (below  approximately  100)  or  if  all  the  included  variable  coef- 
ficients have  F-values  below  100,*  that  equation  also  can  reasonably  be  eliminated 
in  the  maximum  likelihood  optimization  procedure. 

Further  development  of  the  parameter  subset  for  maximum  likelihood  esti- 
mation based  on  regression  estimate  F-values  is  also  desirable.  Experimentation 
reveals  that  both  less  computation  and  better  estimates  result  if  maximum  likeli- 
hood is  used  to  refine  estimates  of  only  those  parameters  whose  F-value  is  more 
than  4%  of  the  maximum  parameter  F-value  for  a particular  equation.  The  4% 
cutoff  corresponds  to  parameter  estimate  confidence  five  times  less  than  those 
of  the  "best"  parameter  estimate  (c.f.  Table  t . 3)  . Attempts  to  identify  more 
parameters  with  the  maximum  likelihood  procedure  may  over -parameterize  the 
fit  to  the  given  data,  resulting  in  poorer  overall  estimates  of  the  parameters  when 
compared  to  the  estimates  obtained  with  a smaller  parameter  set. 

If  the  parameters  remaining  after  insignificant  modes  have  been  eliminated, 
the  regression  subset  is  decomposed  into  two  parts.  The  first  part  is  those  terms 

An  F-value  of  100  corresponds  to  a 95%  parameter  estimate  ccfhfidence  (2  sigma) 
of  about  +20%. 
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which  are  "significant"  by  the  various  F-tests.  The  regression  estimated  values 
of  these  terms  are  defined  more  precisely  by  the  maximum  likelihood  algorithm. 

The  second  part  consists  of  terms  which  are  "insignificant" , but  for  which  numerical 
values  have  been  estimated.  Because  these  latter  estimates  do  constitute  additional 
information  about  the  system,  they  are  used  as  fixed  values  in  the  maximum  likeli- 
hood algorithm.  Hence,  all  regression  estimates  are  used  for  the  maximum  likelihood 
identification;  one  subset  as  a priori  start-up  values  and  the  other  subset  as 
fixed  values.  The  reason  for  including  the  non-zero  values  of  the  "insignificant" 
subset  is  that  these  regression  estimates,  though  not  perfect,  should  be  closer 
to  the  actual  value  than  an  estimate  of  zero.  The  estimates  of  "identified"  parameters 
improve  if  the  estimates  of  "non-identified"  parameters  are  better. 

Overall,  the  intent  of  emphasizing  reduction  of  the  parameter  set  for  maximum 
likelihood  estimation  is  twofold.  First,  and  most  importantly,  the  problem  of 
over-parameterization  and  the  attendant  risk  of  poorer  general  quality  of  all 
parameter  estimates.  This  difficulty  arises  largely  because  with  a finite  set 
of  noisy  data,  the  addition  of  more  parameters  to  obtain  a better  model  fit  to  the 
data  will  eventually  result  in  the  model  attempting  to  fit  the  noise  process  at 
the  expense  of  knowledge  of  the  actual  underlying  process  of  interest.  Thus  , 
the  attempt  to  identify  large  numbers  of  parameters  merely  for  the  sake  of  the 
accomplishment  may  be  counter-productive.  The  second  reason  is  that  of  efficiency. 
Diminishing  returns  result  in  terms  of  estimate  confidences  by  increasing  the 
parameter  set  identified  from  a fixed  data  set,  but  the  computation  effort  required 
continues  to  increase  at  least  as  fast  as  the  parameter  set  size. 

Table  6.4  presents  the  parameter  sets  for  the  two  demonstration  examples, 

lateral  case  1 and  longitudinal  case  2.  Included  in  the  sets  used  for  maximum 

likelihood  estimation  in  addition  to  the  parameters  selected  by  statistical  means 

are  the  constant  terms  for  each  state  equation  and  the  bias  on  sideslip  angle,  p, 

measurement.  These  additional  parameters  do  not  have  the  F-value  statistic 

for  significance  ranking,  but  are  nevertheless  felt  to  be  potentially  important 

enough  to  be  included.  In  the  case  of  the  constant  terms,  the  value  is  that  of 

the  primary  aerodynamic  coefficient  values  (C„  ,C  ,C  , C , C , C ) at  the 

x m w x y z 

mean  a and  p for  the  experiment  data.  Thus,  the  final  rule  in  determining  the 
parameter  set  should  always  be  the  judgment  of  the  user. 
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TABLE  6 . 4 

Parameter  Set  Decomposition 
Lateral  Response  Case  1 


Parameter  Estimates 
Refined  by  Maximum  Likelihood 


Other  Non-Zero  Regression 
Parameter  Estimates  Used  in 
Maximum  Likelihood 


£ 

n 

a 

ap 

£ 2n 
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Longitudinal  Response  Case  1 


Parameter  Estimates 
Refined  by  Maximum  Likelihood 


Other  Non-Zero  Regression 
Parameter  Estimates  Used  in 
Maximum  Likelihood 


Jm  2 2 
a p 


'’m  R 2 

8sa 
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6. 2. 3. 2 Evaluating  Results  of  Maximum  Likelihood  Estimation 

In  addition  to  the  actual  parameter  estimates  obtained  using  the  maximum 
likelihood  procedure,  covariances  oi  the  estimates  are  available  from  the  information 
matrix  (see  Appendix  A)  using  the  sensitivities  (of  the  other  variables)  in  the 
model  to  the  parameters  gg,  ggj,  covariances  as  well  as  time  histories 

can  also  be  propagated  for  estimates  of  states,  measurements,  and  nonlinear 
aerodynamic  coefficients.  When  plotted  with  actual  time  histories  available  from 
the  simulated  data,  the  estimated  time  histories  and  associated  covariances  form 
an  excellent  method  for  evaluating  the  results  of  the  entire  identification  procedure. 

Plots  of  only  the  actual  and  estimated  time  histories  are  useful,  but  the 
addition  of  covariances  for  the  estimates  is  helpful  in  examining  poor  time  history 
fits.  Especially  in  the  case  of  aerodynamic  coefficients  composed  of  several 
parameters,  high  covariances  resulting  from  low  information  content  in  the  data 
should  still  bracket  most  of  the  actual  variation  in  a coefficient,  even  if  the  estimated 
value  does  not  conform  very  well. 

Some  latitude  must  be  allowed  in  the  covariance  time  histories,  however, 
since  the  available  covariances  are  also  only  estimates  of  the  actual  covariances . 

The  estimated  covariances  will,  in  general,  be  too  small  for  several  reasons: 

1.  Covariances  are  generated  based  on  the  information  matrix  for  only 
the  parameters  in  the  maximum  likelihood  search.  Inclusion  of 
less  significant  parameters  in  the  information  matrix  calculation 
would  tend  to  increase  the  covariances  of  all  parameters.  The  amount 
of  the  increase  depends  on  significance  of  the  parameters  added. 

2.  Uncertainty  in  other  fixed  constants  not  estimated  in  the  model  such 
as  mass,  inertias,  engine  thrusts,  and  instrument  misalignments 
are  not  accounted  for  in  the  information  matrix.  Inaccuracies  in- 
herent in  these  estimates  as  well  as  uncertainty  in  the  states  incor- 
porated as  time  histories  will  also  increase  covariance  estimates. 
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3.  Finally,  estimated  covariances  are  linearized  estimates  which  pre- 
sumes all  model  parameters  are  at  their  "actual"  values.  In  fact, 
the  modeling  process  itself  involves  simplifications  which  make 
it  impossible  for  the  model  parameters  to  always  equal  the  "real" 
values.  These  modeling  errors  are  manifested  both  as  biased  parameter 
estimates  and  conservatively  small  covariance  estimates. 

For  the  output  error  form  of  the  maximum  likelihood  estimation  procedure, 

2 

linearized  covariance  estimates  for  the  aerodynamic  coefficients,  o^,  are, 


2 _ dc  . .-1  dc 

Oc  - diag  j-0M  dl- 


for  the  state  estimates  , o , 

-x 


Jl  dx  ..-1  dx 

a = diag  -tq 

-x  b dO  - d9 


(6.4) 


(6.5) 


and  the  measurement  estimates, 


4 = diag 


(6.6) 


In  Eq.  (6.6),  the  covariances  due  to  estimates  of  measurement  noise,  implicit  in 
the  R matrix,  are  incorporated  additivcly  to  the  covariance  due  to  the  other  para- 
meter estimates.  This  formulation  is  a result  of  assuming  the  elements  of  R 
are  uncorrelated  with  other  parameter  estimates.  These  simplifying  assumptions 
considerably  reduce  the  computational  demands  of  the  maximum  likelihood  function 
optimization  while  leaving  unaffected  the  final  parameter  estimates.  Not  making 
this  assumption  will  result  in  lower  covariance  estimates  on  predicted  outputs. 
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6 . 2 . 3 . 3 Demonstration  Examples 


Table  6 . 5 lists  a priori  maximum  likelihood  parameter  values  (obtained  from 
the  regression  algorithm  results),  final  parameter  estimates,  parameter  F-values, 
and  other  optimization  results  for  the  two  example  cases.  The  F-value  presented 
for  maximum  likelihood  parameter  estimates  is  calculated  in  the  same  way  as  it  is 
for  the  regression: 


F 


(6.7) 


The  F-value  associated  with  the  a priori  parameter  estimate  is  given  for  both  the 
regression  and  zero^  iteration  of  the  maximum  likelihood  for  comparison  purposes , 
Comparison  of  these  tabular  results  reveals  the  following  important  points: 

a.  Final  maximum  likelihood  F-values  are  higher  than  a priori  values, 
indicating  improvement  in  estimate  quality  over  regression  results. 


b.  F-value  ordering  indicates  relative  parameter  significance  remains 
close  to  regression,  lending  confidence  to  using  the  regression  for 
parameter  set  specification. 

c.  F-values  from  the  regression  and  a priori  maximum  likelihood  results 
are  of  the  same  order  of  magnitude. 

Figures  6.3  through  6.8  are  time  history  plots  of  measurements  and  various 
aerodynamic  coefficients  showing  actual  values  and  estimated  values  with  +2o 
confidence  limits.  The  most  prominent  characteristic  of  both  demonstration  cases 
is  the  excellent  measurement  estimate  fits  to  the  actual  data.  Estimate  2o  confidence' 
limits  bracket  the  measurement  noise  induced  variations  of  the  actual  data. 

Estimate  fits  to  the  nonlinear  aerodynamic  coefficients  reveal  a range  from 
excellent  estimates  to  fairly  significant  biased  estimates.  For  the  lateral  motion 
case,  the  major  coefficient  estimates,  , and  are  very  good;  again, 

with  confidence  limits  bracketing  actual  data  most  of  the  time.”  Estimates  of  control 
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TABLE  6.5 

Maximum  Likelihood  Results  for  Demonstration  Cases 

\ 

Lateral  Response  Case  1 

J = 15.036 

Measurement  Noise  Covariance 

o 

J, = 18.543 

Estimates: 

f 

Jf'Jf-l/Jf  = 2-125  x 

10"7 

P 

3.38  x 10  radians 

-3 

No . of  Iterations  : 

= 5 

y 

9.40  x 10  g 
_ 0 

p 

7.05  x 10  radians/sec 

r 

-3 

5.32  x 10  radians/sec 

A PRIORI 

FINAL 

PARA- 

ESTIMATE 

F -VALUE 

1 

METER 

(FROM 

REGRESSION) 

REGRESSION 

MAX.  LIKELIHOOD 

ESTIMATE 

F-VALUE 

CS 

P 

-6.97  x Hf2 

3.24  x 103 

7.10  x 103 

-6.73  x 10'2 

6.77  x 104 

C 

-4.38  x 10~2 

3.36  x 303 

3.26  x 102 

-4.95  x 10'2 

5.6]  x 103 

n6 

r 

ct 

6 

-3.01  x 10'2 

3.29  x 102 

3.57  x 103 

-2.68  x 10'2 

3.12  x 103 

a 

ct 

2.02  x 10'3 

— 

3.65  x 102 

1.42  x 10'3 

1.85  x 103 

o 

2.45  x I01 

2.25  x 101 

c« 

ra 

1.49  x I02 

2.92  x 102 

1.76  x 103 

C 

-5.06  x 10_1 

1.32  x )03 

1.28  x 103 

-5.36  x I0"1 

1.54  x 103 

VP 

C 

4.64  x 10-2 

1.56  x 102 

5.67  x 105 

5.88  x 10'2 

4.01  x 102 

P 

bp 

0.0 

— 

0.0 

6.R9  x I0-3 

1.66  x 102 

c 

] ?A  v in'2 

7 14  v in1 

2.07  x I01 

6.87  x 10'3 

9.38  x 101 

n6 

a 

Cn 

1 o 

-4.04  x 10‘4 

— 

5.67 

3.29  x 10‘4 

7.65  x 101 

6.20  x 101 

2.7  x 101 

1.55  x 101 

cnp3 

2.47 

1.09 

p 

c 

3.69  x )0'4 



9.36  x 10'2 

-3.40  x 10"3 

1.16  x 101 

y° 

i m 


TABLE  6.5  (Continued) 

Maximum  Likelihood  Results  for  Demonstration  Cases 
Longitudinal  Response  Case  1 


J = 14.965 
o 

Jf  = 18.397 

Jf~Jf-l/Jf  = 5,476  x 10" 
No.  of  Iterations  = 7 


Measurement  Noise  Covariance 
Estimates: 

-3 

a 1.94  x 10  radians 

x 6.05  x 10  5 g 

z 1. 01  x 10~2  g 

_2 

q 1.17x10  radians/sec 


PARA- 

METER 


ESTIMATE 

(FROM 

REGRESSION) 


-7.60  x 10_1 
-2.82 

-4.53  x 10'1 

-4.42  x 10'2 
5.52 
4.57 

-8.40  x 102 
-5.06  x 10*2 
-1.12 

-8.49  x 101 
-1.49  x 10*2 


A PRIORI 


F -VALUE 


REGRESSION  MAX.  LIKELIHOOD 


1.71  x 10 


8.50  x 10 


3.59  x 10‘ 


2.75  x I0J 


1.58  x 10‘ 


5.20  x 10 


<;  i U,  v in‘ 


1.02  x 10‘ 


8.60  x 101 


3.10  x 10J 
7.26  x 102 

1.01  x 103 

7.30  x 102 

2.01  x 102 
2.53  x 101 

2.89 

6.45 

1 . 2?  x 3 o2 
1.62 

4.23  x 101 


FINAL 

ESTLMATE 

F-VALUE 

-7.62  x 10'1 

3.27  x 105 

-2.86 

1.19  x )05 

-4.64  x 10’1 

1.76  x 104 

-4.04  x 10~2 

9.50  x I03 

-5.79 

6.87  x 103 

4.45 

2.89  x 103 

-6.11  x 102 

8.61  x 102 

-7.80  x I0'2 

8.45  x 102 

5.4*.  x !02 

-6.73  x 101 

3.99  x 102 

-1.43  x 10'2 

1.42  x 102 

8.86  x 10-2 

9.00  x 101 

I 


NONE 


Pftf 


Figure  6 .5  Lateral  Control  Effectiveness  Coefficient  Response  (Simulated  or  Estimated)  to  Combined 
Aileron/Rudder  Input  (Lateral  Response  Case  1,  Ref.  Fig.  6.2a). 


<9>N0J1W831300U  TtfO!l83A 


<9>N01JU83*133DW  TW  l OfUl  9N03 


Longitudinal  Measurement  Response  (Simulated  and  Estimated)  to  Sinusoidal  Stabilator  Input 
(Longitudinal  Case  1,  Ref.  Fig.  6.2c). 


Longitudinal  Static  Coefficient  Response  (Simulated  and  Estimated)  to  Sinusoidal  Stabilator 
£nput  (Longitudinal  Case  1,  Ref.  Fig.  6.2c). 


Response 

Estimated  (with 
2<7  bounds) 


derivatives  such  as  C , C and  C0  tend  to  be  biased,  however.  The  two 

n~  n~  Xo 

or  oa  oa 


prime  reasons  for  the  biased  estimates  are: 


Modeling  errors  due  to  higher  order  variation  of  actual  coefficients 
with  a and  P than  allowed  for  in  coefficient  model  polynomial  expansion. 


A small  amount  of  information  about  control  der  ivatives  is  available  as 


a result  of  the  short  time  of  control  application. 


In  the  longitudinal  case,  very  good  estimates  of  C^,  and  are  produced, 
but  Cm  is  in  error,  especially  at  the  extremes  of  its  variation.  Also,  the  biases 
in  the  other  coefficients  are  again  evident. 


The  reason  for  the  degraded  estimates  of  both  C and  C is  primarily 

m m „ 

os 


due  to  lack  of  information  about  the  highly  nonlinear  portion  of  the  vs.  a 
curve  over  which  the  experiment  is  centered  (see  Fig.  (6.9))  . Information  about 
C is  small  with  the  high  level  sinusoidal  variation  of  stabilator  force  throughout 


the  experiment.  The  errors  in  the  C estimates,  in  turn,  detract  from  the  C 

m mt 


estimate  since  an  incorrect  estimate  of  will  be  compensated  for  by  an  altered 


C estimate  so  that  measurement  time  histories  will  still  match.  This  situation 
mc 
os 


requires  data  with  better  information  to  improve  estimates  of,  and  thereby  distinguish 
between,  the  two  coefficients , as  is  shown  in  Section  6.4.  Input  design  is  thus 
seen  to  be  extremely  important  in  generating  the  maximum  information  for  estimating 
all  important  parameters  simultaneously,  since  parameter  estimate  interdependence 
results  in  the  poor  estimate  of  one  parameter  causing  poor  estimates  for  other 
parameters . 


The  basic  features  and  implementation  of  the  overall  parameter  estimation 
procedure  are  demonstrated.  Excellent  reproduction  of  measurement  time  histor- 
ies are  achievable;  and  within  the  limitations  of  discussed  of  modeling  and  avail- 
able information  in  the  data,  linear  and  nonlinear  aerodynamic  coefficients  can  be 
identified  extremely  well.  In  the  following  sections,  further  exploration  of  the  capa- 
bilities of  the  identification  procedure  are  demonstrated.  Included  are  the  effects 
of  measurement  noise  levels,  process  noise,  input  variation,  and  data  length. 


Finally,  a comprehensive,  full  six  degrees-of-freedom , 31  parameters , 15  second 
data  length  case  is  presented.  In  the  truest  test  of  an  estimation  procedure,  the 
estimated  parameters  from  this  final  case  are  used  to  successfully  predict  time 
histories  for  entirely  different  control  inputs. 

6 . 3 NOISE  EFFECTS 

In  general,  the  effect  of  increased  noise,  whether  process  or  measurement, 
on  parameter  estimation  is  to  decrease  confidence  in  the  estimates.  This  section 
describes  experiments  to  illustrate  the  impact  of  noise  on  the  integrated  parameter 
identification  procedure.  Results  show  that  even  for  quite  high  levels  of  measure- 
ment and  process  noise,  the  identification  procedure  is  still  viable  and  yields 
reliable  parameter  estimates. 

6.3.1  Measurement  Noise 

The  base  demonstration  cases  of  Section  6.1  employ  simulated  data  in  which 
thei'e  is  instrumentation  measurement  noise  present.  A median  level  for  those 
noises  is  used  based  on  information  in  ref.  [1]  . Comparison  response  data, 
exactly  the  same  as  the  lateral  motion  response  case  1 except  for  higher  measure- 
ment noise  levels,  is  generated  by  simulation  in  order  to  assess  the  effect  of 
measurement  noise.  The  experiment,  lateral  response  case  2,  is  described  in 
Figure  6.2,  and  values  are  given  in  Table  6.6  for  both  base  and  high  measure- 
ment noise  parameters  from  the  instrumentation  simulation  equations  of  Appendix 


Table  6.7  compares  regression  results  for  the  high  measurement  noise 
case  to  those  obtained  previously  for  base  measurement  noise.  As  expected, 
there  is  a substantial  decrease  in  included  parameters  due  to  the  higher  coefficient 
covariances  (lower  F-values)  that  result  from  the  increased  noise.  For  what 
otherwise  is  the  same  data,  the  addition  of  more  measurement  noise  reduces  the 
amount  of  useful  information  for  identifying  parameters. 

Significantly,  however,  the  parameters  that  the  regression  does  identify 
in  the  high  measurment  noise  case  are  exactly  the  same  as  the  most  important 
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TABLE  6.7 

Regression  Results  for  High  Measurement  Noise  Case 


OVERALL  FIT 
F-RATJO 

PARAMETER 

PARAMETER 

ESTIMATE 

Cv 

-4.83  x 10'1 

445.9 

YP 

C*o 

1.20  x 10"3 

% 

-8.11  x 10'2 

c£ 

-3.16  x 10'2 

962.7 

6 

a 

cv 

-1.11  x 103 

o 

K> 

0 

9.41  x 10"4 

S 

r 

-4.11  x 10"2 

% 

4.07  x 10'2 

473.9 

c»„ 

~a 

1.77  x 10*2 

cv 

2.14 

Cn 

o 

1.40  x 10’3 

PARAMETER 

F-VALUE 

4.46  x 102 

2.32  x 103 
1.89  x 102 

5.48  x 101 

3.97  x 102 

8.10  x 101 
5.66  x 101 

2.84  x 101 


I 


parameters  isolated  by  the  regression  for  the  base  measurement  noise  case. 

This,  despite  the  fact  that  there  is  measurement  noise  in  this  formulation  of  the 
regression,  tends  to  bias  parameter  estimates.  The  example  thus  reinforces  confi- 
dence in  the  regression's  primary  purpose  of  isolating  the  most  important  para- 
meters from  the  full  admissible  set,  even  in  the  presence  of  severe  instrumentation 
errors . 

Maximum  likelihood  refined  parameter  estimates  for  the  high  measurement 
noise  case  are  given  in  Table  6.8.  Note  that  the  final  F-values  are  approximately 
an  order  of  magnitude  below  those  for  the  bias  measurement  noise  case.  Also, 
the  final  F-valuc  ordering  is  considerably  changed  from  the  a priori  ordering, 
indicating  severe  estimate  biases  from  the  regression. 
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TABLE  6.8 

Maximum  Likelihood  Results  for  High  Measurement  Noise  Case 

Measurement  Noise  Covariance 
Estimates: 

P 5.24  x 10  ^ radians 

y 1.45  x 10  ^ g 
-2 

p 3.64  x 10  radians/sec 

_2 

r 1.47  x 10  radians/sec 


? 

I 

■ i 

i I 

: 
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Comparison  of  the  parameter  estimates  from  bias  and  high  measurement 

noise  cases  on  a one-for-one  basis  docs  not  indicate  very  close  correspondence 

in  most  cases.  Only  C , b„ , C , and  C for  the  high  measurement  noise 

Vp  P n6a  y 

case  have  estimates  within  2o  of  the  base  measurement  noise  estimates.  However, 
because  the  two  parameter  sets  are  not  completely  analogous,  it  is  only  legitimate 


A PRIORI 

FINAL 

PARA- 

ESTIMATE 

F -VALUE 

METER 

REGRESSION 

MAX.  LIKELIHOOD 

ESTIMATE 

F-VALUE 

2 o CONFIDENCE 

% 

-8.11  x 10*2 

2.32  x 103 

9.01 x 103 

-7.98  x 10*2 

8.03  x It3 

+.18 x 10*2 

Cn 

o 

1.40  x 10"3 

— 

1.53  x 101 

1.97  x 10*3 

2.55  x 103 

+.08 x 10*3 

% 

r 

-4.11  x 10*2 

3.97  x 102 

1.10  x 102 

-6.93  x 10'2 

2.39  x 103 

+.28  x 10'2 

% 

4.07  x 10'2 

8.10  x 101 

2.53 

1.67  x 10'1 

7.11  x 102 

+.13  x 10'1 

% 

-4.83  x 10*1 

4.46  x 102 

2.51 x 102 

-5.07  x 10'1 

5.43  x 102 

+.44  x 10* 1 

V 

2.14 

2.84  x 101 

5.47  x 10*1 

-1.05  x 101 

2.58  x 102 

+.13  x 101 

C*6 

a 

-3.16  x 10*2 

1.89  x 102 

6.07  x 101 

-2.18  x 10*2 

2.03  x 102 

+.31  x 10*2 

bP 

0.0 

— 

0.0 

8.53  x 10*3 

6.60  x 101 

+2.1  x 10*3 

C*o 

-9.41  x 1C*4 

— 

4.16 

1.69  x 10*4 

6.10 

+1 .4  x If**4 

\ 

a 

1.77  x 10'2 

5.66  x 101 

1.03  x 101 

6.11 x 10"3 

5.45 

+5.3  x 10*3 

C'o 

-1/20  x 10'3 

— 

2.22  x 10'1 

-2.56  x 10*3 

2.59 

♦3.2  x 10*3 

J = 11.693 
o 

Jf  = 15.016 

Jf"Jf-l/Jf  = 4-729  x 10 
No.  of  Iterations  = 7 
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to  compare  the  net  nonlinear  aerodynamic  coefficients  composed  of  several  parameters 
since  nearly  similar  fits  can  be  obtained  with  different  parameter  sets.  Figures 
6.10  through  6.12  are  the  time  history  comparisons  for  the  high  measurement 
noise  case,  including  the  aerodynamic  coefficient  traces.  Of  all  the  coefficient 
estimate  time  histories,  only  that  of  is  seriously  different  from  the  actual 
time  history.  Thus,  despite  the  large  variations  in  specific  parameter  estimates 
between  high  and  base  measurement  noise  cases,  the  variation  estimates  of  the 
total  coefficients  are  not  comparatively  in  error.  In  fact,  the  primary  difference 
is  the  expected  increase  in  covariances  for  the  high  noise  case. 


Time  histories  for  measurement  estimates  again  exhibit  excellent  fit  to 
the  actual  data.  These  consistently  good  fits  of  estimated  and  actual  measurements 
are  indicative  of  the  basic  principle  that  for  a sufficiently  generalized  model 
(such  as  the  one  used  here)  and  enough  independent  parameters  (legrees-of- 
freedom)  almost  any  measurement  data  can  be  matched  very  closely.  The  difficulty 
is  that  too  close  a measurement  match  will  result  in  decreased  relevence  of  the 
parameter  values  to  the  characteristic  they  are  supposed  to  model.  This  can 
be  seen  in  the  yaw  rate  time  history  and  the  coefficient.  Close  inspection 
of  the  yaw  rate  estimate  reveals  tracking  of  the  actual  yaw  rate  value  to  the  point 
of  following  some  of  the  noise  variations.  This  effect  has  been  discussed  in  Section 
3.3.2  where  it  is  noted  that  one  effect  of  choosing  more  parameters  than  actually 
required  is  to  allow  the  estimates  to  ''track  the  noise" . For  this  data,  with  the 
high  noise  present  masking  to  a large  extent  higher  order  variations  in  C^,  the 
attempt  to  estimate  both  C , and  C is  probably  over-parameterization, 

np3  np 

The  parameter  set  specification  described  in  Section  6. 1.2  eliminates  a 
great  deal  of  the  risk  of  over-parameterization.  In  some  cases,  the  regression 
technique  will  tend  to  select  higher  order  polynomial  terms  when  lower  order 
terms  would  be  satisfactory,  especially  when  high  noise  levels  are  present. 

The  regression  technique  has  an  option  for  forcing  in  lower  order  terms  which 
is  used  in  such  cases.  Such  "forcing  in"  does  not  necessarily  require  that  the 
terms  be  included,  only  that  they  not  be  deleted  unless  found  highly  insignificant 
relative  to  other  terms. 
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Figure  6.11  Lateral  Static  Coefficient  Response  (Simulated  and  Estimated)  for  Measurements  with  Increased 
Instrument  Noise  (Lateral  Case  2,  Ref.  Fig.  6.2a)  . 


Lateral  Control  Effectiveness  Coefficient  Response  (Simulated  and  Estimated)  for  Measurements 
with  Increased  Instrument  Noise  (Lateral  Case  2,  Ref.  Fig.  6.2a)  . 


6.3.2  Process  Noise 


The  formulation  of  the  maximum  likelihood  estimation  procedure  used  is 
the  output  error  form  which  does  not  identify  the  process  noise  in  the  experiment 
data.  The  chief  advantage  of  this  formulation  is  a large  reduction  in  the  computations 
needed  to  carry  out  the  nonlinear  estimation  process.  Unless  gusts  are  specifically 
to  be  estimated,  such  computational  increases  may  not  be  necessary.  In  addition, 
properties  of  the  regression  program  reduce  the  effect  of  such  noise  on  the  estimate. 

Even  when  no  process  noise  is  assumed,  some  process  noise  type  effects 
always  exist  as  a result  of  modeling  errors.  In  order  to  evaluate  the  identification 
procedure's  robustness  to  a more  significant  level  of  process  noise,  an  experiment 
is  simulated  in  which  moderate  levels  of  white  wind  gust  turbulence  are  added 
to  ' ertical  and  horizontal  aircraft  velocities  for  the  lateral  mode  demonstration 
case  (base  measurement  noise  is  also  present)  . Figures  6.13  through  6. 14  show 
the  measurement  and  coefficient  time  histories.  Table  6.8  gives  the  process 
noise  parameter  values,  and  the  experiment  is  described  under  lateral  response 
case  3 in  Figure  6.2. 


Regression  results  are  given  in  Table  6.9  and  maximum  likelihood  esti- 
mates in  Table  6.10.  C , C«  , C , and  C are  once  again  isolated  by  the 

y P n8r  np 

regression  as  the  most  significant  parameters  for  the  three  equation  fits  (as  in 

the  base  and  high  measurement  noise  cases)  . Beyond  those  four,  only  C , 

P 

corresponds  with  significant  sets  for  those  other  cases . That  leaves  five  different 

parameters  picked  in  the  process  noise  case.  The  aerodynamic  coefficient 

estimates  appear,  as  well  as  can  be  determined  with  the  high  level  of  noise  in 

the  actual  data,  excellent,  with  the  possible  exception  of  C . The  2 a covariance 

8r 

estimates  on  the  aerodynamic  coefficients  are,  however,  considerably  less  than 
the  apparent  noise  in  the  actual  data.  This  is  because  the  maximum  likelihood  out- 
put error's  structure  does  not  account  for  the  process  noise  responsible  for  the 
random  variation  in  coefficients.  Instead,  the  process  noise  is  similar  in  effect 
to  additional  measurement  noise,  as  can  be  seen  in  the  measurement  estimate 
covariances  which  bracket  the  acUial  measurement  variations  due  to  both  measure- 
ment and  process  noise.  The  increased  uncertainty  in  the  net  measurements 
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TABLE  6.9 

Regression  Results  for  Process  Noise  Case 


EQUATION 


OVERALL  FIT 
F-RATIO 


PARAMETER 


PARAMETER  PARAMETER 
ESTIMATE  F-VALUE 


-5.25  xio'1  9.68  x 102 
-6.43  x 10'3  


4.07  x 10  -1  6.26  x KT 


9.83  x 10  * 8.39  x 1(T 


4.60  x 10 


-3.71  x 10 


1.23  x 10  J 2.61  x KT 


3.66  x 10 


is  reflected  in  aerodynamic  coefficient  covariances  higher  than  for  both  of  the 
measurement  noise  cases.  The  covariances  are  not  as  large,  however,  as  would 
be  the  case  if  proce  ss  noise  was  taken  directly  into  account.  Nevertheless,  the 
actual  coefficient  estimates  themselves  are  still  very  good. 

6.4  METHODS  FOR  IMPROVING  PARAMETER  ESTIMATES 


The  ultimate  limit  of  the  accuracy  of  any  parameter  estimate,  compared 
to  its  "true"  value,  is  set  by  the  relative  effect  the  parameter  has  on  the  response 
data.  This  is  determined  by  the  overall  information  content  of  the  data.  In  general, 


c 


TABLE  6.10 

Maximum  Likelihood  Results  for  Process  Noise  Case 


J = 10.659 
o 

Jf  = 12.028 

Jf“Jf-l/Jf  = 7-657  x 10' 
No.  of  Iterations  = 7 


"Measurement"  Noise  Covariance 
Estimates 

P 5 . ?f  x 10  2 radians 

y 2.92  xlO"2 

-2 

p 9.49  x 10  radians/sec 

-3 

r 5.64  x 10  radians/sec 


A PRIORI 

FINAL 

PARA- 

r-v/ 

vLL  2 

METER 

ESTIMATE 

REGRESSION 

MAX . LIKELIHOOD 

ESTIMATE 

F-VALUE 

2o  CONFIDENCE 

C£ 

P 

-1.02  x 10"1 

6.26  x 102 

9.68  x 102 

-8.60  x 10'2 

1.66  x 103 

+.42  x 10-2 

C 

n6 

r 

-3.71  x 10'2 

4.60  x 102 

1.41  x 102 

-4.99  x 10‘2 

9.45  x 102 

+.32  x 10'2 

c 

n 

r 

-4.26  x 10'1 

4.35  x I01 

1.02  x 101 

-1.26 

2.25  x 102 

+ .17 

C 

n 

-5.57  x 10_1 

1.29  x 102 

6.92 

-1.50 

1.50  x 102 

+ .24 

°p 

c 

np 

4.66  x 10  ^ 

2.26  x 102 

8.89 

1.25  x 10"1 

1.01  x 102 

+.24  x 10'1 

9.00  x 10_1 

8.39  x 101 

1.47  x 102 

6.71  x 10*1 

8195  x 101 

+1.4  x 10'1 

°p 

% 

-5.25  x 10'1 

9.68  x 102 

4.98  x 101 

-5.59  x 10'1 

3.95  x 101 

+1.8  x 10"1 

C8 

4.60  x 10~4 

— 

1.68 

9.22  x 10"4 

1.91  x 101 

+4.2  x 10'4 

O 

c 

1 .23  v 10_1 

2.61  v 101 

1.67  x 10*3 

-9.56  x ]0-1 

1.58  v 101 

+4.8  x 10_1 

c , 

np3 

2.10 

4.90  x 101 

2.16  x 10*1 

-2.72  x 101 

1.45  x 101 

+1.4  x 101 

C. 

n6 

r« 

2.63  x 10_1 

3.12  x 101 

5.18 

1.48  x 10'1 

5.34 

+1.2  x 10'1 

c 

n 

3.66  x 10~4  ' 

— 

7.84  x 10'1 

5.00  x 10'4 

3.55 

+5.3  x 10'4 

o 

v 

-6.43  x 10'3 
0.0 

2.75 

0.0 

-4.53  x 10‘3 
-4.62  x 10'3 

1.83 

6.32  x 10*1 

+6.7  x 10'3 
+12.0  x 10'3 

parameter  estimates  can  be  only  as  "good"  as  this  data  information  content. 

As  shown  in  Section  6.3,  increasing  noise  in  the  data  reduces  the  quantity  of 
useful  information  available  for  parameter  identification.  This  section  discusses 
practical  means  for  improving  the  useful  information  and  thus  providing  more 
accurate  estimates. 


The  following  methods  are  considered  for  this  improved  jdentifiability 
objective: 

a.  Input  design 

b.  Data  aggregation  (use  of  a priori  information) 

c.  Data  length  increase  (use  of  more  data) 

d.  Addition  of  more  parameters 

6.4.1  Input  Design 


In  the  discussion  of  the  longitudinal  motion  demonstration  example  results, 
the  importance  of  the  input  is  pointed  out  with  regard  to  identifying  all  significant 
parameters  well.  If,  as  in  the  case  with  C . one  parameter  is  poorly  estimated 
as  a result  of  information  unavailable  due  to  an  inappropriate  input,  other  para- 
meter estimates,  C , may  also  be  adversely  affected. 

8 s 

In  an  attempt  to  correct  the  deficiencies  of  that  first  longitudinal  motion 

case,  another  case,  identical  in  all  respects  except  for  the  input,  is  performed. 

For  this  example,  detailed  under  longitudinal  case  2 in  Figure  6.2,  the  stabilator 

input  amplitude  is  reduced  to  2°  in  an  attempt  to  allow  the  influence  of  the  C 
. . m 

nonlinearity  to  have  more  impact  on  aircraft  motion  and  thereby  increase  information 

in  the  data  about  C . 

m 

Results  of  this  case  are  contained  in  Tables  6.11  and  6.12,  and  Figures 

6.15  through  6.18.  For  this  modified  input,  the  C estimate  time  history  has 

a much  better  fit  to  the  variations  in  the  actual  curve,  but  there  is  some  bias. 

The  bias  is  caused  by  the  biased  estimate  of  C apparent  in  Figure  6.16. 

8s 

The  stabilator  input  selected  for  this  case  is  too  low  in  amplitude  to  satisfactorily 

estimate  C . Aircraft  motion  is  primarily  affected  by  the  nonlinearity  in  C 
8s  m 
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TABLE  6 . 12 

Maximum  Likelihood  Results  for  Longitudinal  Case  2 


J =-9.7847 
o 

Jf  = -20.556 
Jf-Jf_ l/jrf  = 2.915  x 10 
No.  of  Iterations  = 7 


Measurement  Noise  Covariance 
Estimates: 

-3 

a 1.21  x 10  radians 

x 6 . 30  x 10  3 g 
l 7.10  x 10~3  g 

_3 

q 2.94  x 10  radians/sec 


A PRIORI 

FINAL 

Para- 

meter 

ESTIMATE 

F-VALUE 

REGRESSION 

MAX.  LIKELIHOOD 

ESTLMATE 

F-VALUE 

c 

z 

-9.65  x 10'1 

— 

8.10  x 104 

-9.71  x 10'1 

7.10  x 105 

c . - 

ma3 

-3.33  x 101 

1.29  x 103 

3.98  x 102 

-3.37  x 101 

1.56  x 104 

cz 

-1.51 

1.02  x 104 

1.32  x 103 

-1.57 

1.52  x 104 

a 

c 

m6 

s 

-5.00  x 10_1 

4.19  x 103 

6.88  x 102 

-5.29  x 10'1 

1.12  x 104 

c 

m 

q 

c 

m 

o 

-6.19 

1.91  x 103 

9.50  x 102 

-7.11 

8.18  x 103 

4.70  x 10'2 

— 

5.44  x 102 

-4.47  x 10*2 

7.20  x 103 

C ? 

m2 

-2.41 

4.35  x 103 

1.11  x 102 

-2.49 

6.25  x 103 

r 

3.41  x 10'2 

3.89  x 103 

-3.43  x 10'2 

4.97  x 103 

o 

1 

1.03  x 101 

c , 

2.20  x 103 

2.20  x 30Z 

1.02  x 101 

1.96  x 103 

a 

c 

1.27 

2.54  x 102 

1.77  x 102 

1.21 

2.58  x 103 

\ 

S 

Figure  6.19  compares  estimates  of  C^  for  the  longitudinal  cases  1 and 
2,  In  case  1,  the  fit  is  poor,  but  there  is  less  bias  than  is  present  in  the  better 
nonlinear  convolution  fit  of  case  2.  If  the  information  available  from  both  experi- 
ments could  be  combined,  a much  better  net  fit  would  result. 
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Longitudinal  Measured  and  Estimated  Response  to  Reduced  Amplitude  Sinusoidal  Stabilator 
Input  (Longitudinal  Case  2,  Ref.  Fig.  6.2b). 


Longitudinal  Static  Coefficient  Response  (Simulated  and  Estimated)  to  Reduced  Amplitude 
Sinusoidal  Stabilator  Input  (Longitudinal  Case  2,  Ref.  Fig.  6.2b)  . 


Longitudinal  Dynamic  Coefficient  Response  (Simulated  and  Estimated)  to  Reduced  Amplitude 
Sinusoidal  Stabilator  Input  (Longitudinal  Case  2,  Ref.  Fig.  6.2b)  . 


Longitudinal  Control  Effectiveness  Response  (Simulated  and  Estimated)  to  Reduced  Amplitude 
Sinusoidal  Stabilator  Input  (Longitudinal  Case  2,  Ref.  Fig.  6.2b). 


6.4.2  Information  Aggregation 


Combining  parameter  estimates  to  aggregate  information  from  two  experi- 
ments is  a fairly  straightforward  procedure  which  involves  both  parameter  esti- 
mates and  information  matrices  (see  Chapter  IV)  . A difficulty  can  arise,  however, 
when  nonlinear  aerodynamic  coefficient  expansions  are  identified.  The  expansions 
are  really  only  valid  for  the  operating  range  of  the  experiment  from  which  they 
are  estimated.  Validity  of  the  expansion  fit  outside  of  that  operating  range  is 
very  questionable  since  significant  and  unpredictable  curvature  variations  in 
the  actual  coefficient  function  are  likely. 

It  is  recommended  at  present  that  aerodynamic  coefficient  estimates  be 
combined  from  various  experiments  of  different  operating  ranges  in  order  to 
develop  an  overall  coefficient  function  estimate.  Coefficient  estimates  should 
be  combined,  on  a point-by-point  basis,  only  in  those  regions  where  estimate 
operating  ranges  overlap. 

Ultimately,  it  would  be  desirable  to  accumulate  previous  estimate  information 
and  use  it  in  the  form  of  a priori  information  matrix  values  for  each  successive 
maximum  likelihood  processing  of  new  data.  This  would  result  in  continuing 
improvement  of  the  estimates  since  each  new  set  of  estimates  would  automatically 
incorporate  information  from  all  previous  experiments.  However,  a major  difficulty 
of  implementing  this  procedure  for  the  estimation  of  nonlinear  aerodynamic  coefficients 
lies  in  the  problem  of  combining  previous  parameter  estimates  which  are  only 
valid  over  pari  of  the  operating  range  of  a new  experiment.  One  possible  technique 
is  to  weight  the  information  matrix  as  a function  of  the  distribution  of  the  measure- 
ments for  the  experiment.  This  is  equivalent  to  developing  a chi-square  distribution 
of  variances  for  various  measurement  samples.  The  value  of  the  information 
matrix  would  then  go  to  zero  at  the  limits  of  measurement  excursions,  reflecting 
the  complete  uncertainty  about  extrapolated  aerodynamic  coefficient  values. 

As  an  example  of  the  present  technique  of  a posteriori  combination  of  aero- 
dynamic coefficient  estimates  from  different  experiments  over  different  operating 

DC 

ranges , lateral  response  case  4 results  for  are  combined  with  those  of 

lateral  response  case  1.  Tables  6.13,  6.14  and  Figures  6.20  and  6.21  describe 


138 


m^yrnm 


TABLE  6.14 

Maximum  Likelihood  Results  for  Lateral  Motion  Case  4 


J =10,401 
o 

Jf  = 15.114 

Jf~Jf-l/Jf  = 2-83  x 10_< 
No.  of  Iterations  = 10 


Measurement  Noise  Variance  Estimates: 

p 9.76  x 10  ^ radians 

y 8 . 34  x 10  ^ g 

-2 

p 4.28  x 10  radians/sec 

-2 

r 1.05  x 10  radians /sec 


A PRIORI 

FINA 

L 

PARA- 

F -VALUE 

METER 

ESTIMATE 

REGRESSION 

MAX.  LIKELIHOOD 

ESTIMATE 

F -VALUE 

C 

n6 

r 

-3.60  x 10'2 

9.50  x I02 

1.18  x 103 

-5.35  x 10"2 

1.69  x 103 

c 

yp 

ci 

p 

-5.01  x 10'1 

1.91 x 103 

1.26  x 103 

-5.12  x 10'1 

1.41  x 103 

-7.35  x 10* 1 

1.37  x 103 

2.67  x 102 

-7.58  x 10‘2 

4.11  x 102 

c* 

2.56 

6.49  x 102 

3.79  x 102 

2.00 

2.64  x 102 

°p 

3.65  x 101 

Cn„5 

5.03  x 101 

1.39  x 102 

8.86  x 101 

2.43  x 102 

P 

C« 

o 

8.99  x 10* 1 

1.77  x 102 

5.27  x )0'4 

7.62  x 101 

c 

-1.01  x 10*1 

1.61  x 102 

5.70  x I01 

-6.85  x 10‘2 

3.98  x 101 

p 

3.79  x 101 

3.28  x 101 

CK„3 

2.58 

1.80  x 102 

1.61 

K 

c 

-2.20 

3.97  x 102 

1.99  x 101 

-1.35 

2.92  x 101 

% 

2.02  x 101 

c 

1.43  x 10'  4 

■ 

1.12 

9.86  x 10"4 

O 

c 

na 

C 

2.33  x 10*2 

1.58  x I02 

3.42  x 101 

9.12  x 10‘3 

6.51  x 101 

1.21  x 10*3 

.... 

1.29 

2.14  x 10"3 

4.54 

lateral  response  case  4 results.  Figure  6.22  shows  the  two  estimates  of  ~ < 
the  actual  function  and  the  combined  estimate.  The  combined  estimate  is  abetter 
fit  to  the  actual  curve  in  the  overlapping  region  than  either  estimate  separately. 


Lateral  Measured  and  Estimated  Response  to  a Combined  Longitudinal/Lateral  Input  (Lateral 


individual  Max . Likelihood  Estimates 


6.4.3  Increasing  the  Number  of  Parameters 

A lack  of  sufficient  parameters  is  frequently  a cause  of  significant  errors 
in  an  identification  model  time  history,  compared  to  the  original  system  output. 
Thus,  it  is  often  erroneously  concluded  that  increasing  the  number  of  parameters 
will  improve  the  time  history  match  and  perhaps  even  improve  the  parameter 
estimates . 

Over-parameterization  is  an  attempt  to  derive  more  model  parameter  esti- 
mates from  data  than  are  actually  identifiable  from  the  useful  information  content 
of  the  data.  In  essence,  the  problem  with  over-parameterization  is  that  the 
degree-of-freedom  allowed  the  estimate  or  fit  by  adding  parameters  becomes 
too  large  relative  to  the  quantity  of  the  data.  The  result  is  that  although  there 
ma}  be  a monotonically  improving  fit  to  the  actual  data,  the  statistical  significance 
of  the  fit  improvement  will  eventually  decrease  (this  situation  is  described  by 
the  F-ratio  test  in  the  regression  program)  . 

If  the  inclination  to  add  parameters  to  improve  the  estimate  fit  to  real 
data  is  continued  beyond  the  point  of  maximal  statistical  significance,  the  overall 
confidence  in  estimates  of  parameters  will  also  decline.  A qualitative  reason 
for  this  phenomenon  is  that  the  estimate  model  begins  trying  to  account  for  noise 
in  the  real  data  in  addition  to  modeling  the  underlying  process  of  interest.  In 
effect,  the  result  is  that  a close  to  perfect  fit  of  the  available  data  may  be  obtained, 
but  the  predictive  capability  of  the  over-parameterized  model  is  very  poor. 

Section  6.2  outlines  the  parameter  set  specification  recommended  for  the 
identification  procedure  used  in  order  to  avoid  the  problem  of  over-parameteriza- 
tion. Advantage  is  taken  of  the  readily  available  regression  statistics  related 
to  overall  fit  and  individual  parameter  significances  in  order  to  pick  a realistic 
parameter  set  for  maximum  likelihood  estimation. 

As  an  illustration  of  the  effects  of  an  over-parameterized  maximum  likeli- 
hood estimation,  the  data  from  lateral  case  1 is  reprocessed  with  22  parameters 
identified  instead  of  the  original  12.  Results  are  given  in  Table  6. 14  and  Figures 
6.23  through  6.25.  Time  history  comparisons  with  the  original  12  parameter 
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Figure  6.24  Lateral  Static  Coefficient  Response  (Simulated  and  Estimated  Based  on  22  Parameters)  to  an 
Aileron/Rudder  Input  (Lateral  Case  2,  Ref.  Fig.  6.2a). 


identification  of  lateral  motion  case  1 (figures  6.3  to  6.5)  demonstrate  that  the 

measurement  estimates  are  improved.  For  coefficients,  however,  only  C is 

n6r 

significantly  improved  and  slightly  improved.  Cg  is  about  the  same,  while 

C is  slightly  worse  and  C0  and  C are  much  worse.  The  net  assessment  is 
n6a  oa  y 

then  tli at  the  22  parameter  identification  produces  slightly  poorer  parameter  esti- 
mates than  the  12  parameter  identification.  There  may  be  a better  optimum  set  size 
between  12  and  22  parametcis,  or  even  below  12;  but  the  gain  to  be  had  is  probably 
small  since  the  change  from  12  to  22  is  so  slight.  The  figure  of  12,  however, 
has  the  distinct  advantage  of  having  been  found  by  a single  pass,  repeatable 
procedure,  while  the  unknown  "better"  optimum  set  size  must  bt  found  by  maximum 
likelihood  processing  with  several  iterative  choices  of  set  size. 

(>AA  Example  of  Results  for  an  Increased  Data  Length  and  Coupled 
Mode  Models 

In  general,  identification  of  more  parameters  is  feasible  without  over- 
parameterizing the  problem  only  if  more  information  about  the  parameters  is 
available.  Previous  discussions  of  Section  6.3  describe  the  impact  of  information 
quality  in  a given  quantity  of  data.  In  this  section,  an  example  is  presented 
where  the  data  quality  is  comparable  to  earlier  base  results  (i.e.  , approximately 
the  same  inputs  are  employed  as  in  lateral  response  case  1 and  longitudinal  response 
case  1),  but  the  quantity  of  data  is  increased  because  the  experiment  is  simulated 
for  15.0  seconds  rather  than  10.0  seconds. 


Details  of  this  experiment  are  given  in  Figure  6.2c  under  lateral/longitudinal 
response  case  1.  Both  lateral  and  longitudinal  plane  inputs  are  used  and  para- 
meter estimates  are  obtained  for  all  five  degrees-of-frecdom . Results  are  in  the 
standard  format  of  Table  6.15  and  Figures  6.26  through  6.29. 

There  are  24  parameters  in  the  significant  set  identified  plus  the  7 diagonal 
elements  of  the  R matrix  (which  are  estimates  of  the  measurement  noise  covariances)  , 
giving  a total  of  31  estimates  . The  total  of  estimates  for  the  comparable  longitudinal 
motion  case  1 and  lateral  motion  case  1,  less  the  estimates  related  to  the  u equation, 
is  30.  Tims,  approximately  the  same  number  of  parameters  arc  singled  out  in 
the  combined  case  to  describe  responses  similar  to  the  decoupled  cases.  Examina- 
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Longitudinal  Measured  and  Estimated  Response  with  Inclusion  of  Cross-Coupling  of  Equations 
and  Increased  Data  Length  (Longitudinal/Lateral  Case  1,  Ref.  Fig.  6.2c). 


Lateral  Measured  and  Estimated  Response  with  Inclusion  of  Cross-Coupling  of  Equations  and 
Increased  Data  Length  (Longitudinal/Lateral  Case  1,  Ref.  Fig.  6.2c). 


Longitudinal/Lateral  Static  Coefficient  Response  (Simulated  and  Estimated)  with  Inclusion  of 
tion  Cross-Coupling  and  Increased  Data  Length  (Longitudinal/Lateral  Case  1,  Ref.  Fig.  6.2c 


Lateral  Static  Coefficient  Response  (Simulated  and  Estimated)  with  Inclusion  of  Equation  Cross 
Coupling  and  Increased  Data  Length  (Longitudinal/Late'-al  Case  1,  Ref.  Fig.  6.2c)  . 


TABLE  6.16 

Maximum  Likelihood  Results  for  Lateral/Longitudinal  Motion  Case  1 


J = 25.956 
o 

Jf  = 28.720 

Jf'Jf-l/Jf  = i-09  x 10" 
No.  of  Iterations  = 4 


Measurement  Noise  Covariance  Estimates 

-3 

a 8.06  x 10  radians 

P 3.90  x 10  ^ radians 

y 4 . 52  x 10  k g 

z 3.91xl0_2g 

-2 

p 2.94  x 10  radians/sec 

q 9.32  x 10  ^ radians/sec 

-3 

r 7.46  x 10  radians/sec 


PARAMETER 

A rujoill 

FINAL 

ESTIMATE 

K -VALUE 

ESTIMATE 

F- VALUE 

C 

zo 

-8.24  x 10_) 

1.45  x 104 

-8.06  x IO*1 

3.13  x 104 

C 

mo 

-3.85  x 10‘2 

2.65  x 1C3 

-4.26  x 10~2 

5.09  x 103 

c 

^■6 

s 

-4.24  x 10_I 

2.77  x 103 

-4.60  x 10‘J 

4.98  x 103 

S 

-1.28  x 10'1 

1.33  x 103 

-1.30  x 10'1 

3.91  x 103 

c 

za 

-2.37 

1.15  x 103 

-2.40 

2.67  x 103 

C 

tn 

q 

-5.20 

9.07  x 102 

-6.22 

2.08  x 103 

c 

n6 

r 

-4.77  x 10"2 

8.33  x 102 

. -4.37  x 10~2 

1.75  x 103 

C 

np 

1.08  x 10‘ 1 

6.47  x 102 

8.96  x 10'2 

1.28  x 103 

C8 

P 

-1.59  x 10"1 

3.95  x 102 

-1.76  x JO'1 

1.06  x 103 

% 

-5.51  x 10'1 

6.09  x 102 

-5.50  x 10"1 

9.06  x 102 

C*e 

- a 

-4.51  x IO-2 

4.92  x 102 

-4.01  x 10'2 

8.95  x 102 

-1.90  x 101 

2.48  x 102 

-2.39  x 101 

8.28  x 102 

c 

* o 

-9.60  x IO-3 

4.12  x 102 

-8.69  x 10~3 

5.46  x 102 

c 

n 

o 

1.07  x 10'3 

1.26  x 102 

1.07  x 10'3 

3.61  x 102 

C 

m 

a 

5.58  x 10"2 

2.54  x 101 

1.07 x 10_1 

3.53  x 102 

cc 

r 

-2.77  x 10_1 

1.07  x 102 

-3.10  x 10'' 

3.09  x 102 

ct 

I’u 

9.59  x 10~J 

3.06  x 101 

1.48 

2.53  x iO2 

V 

7.68 

5.21  x 101  • 

5.97 

7.15  x 101 

Cya 

r 

-4.70  x 10'2 

2.65  x 101 

5.99  x 10  2 

6.16  x IO1 

C , 
"«2 

— 

— 

-8.12  x 10'2 

5.21  x 101 

C , 
m 2 

(i 





3.73  x 10’’ 

3.57  x 101 

c, 

O 

2.80  x 10'5 
- 1 

2.67  x 10'1 

1.97  x 10  4 

2.79  x IO1 

c , 

y 

3.56  x 10  1 

2.09  x 10 

) .07  x 10 

5.19 

c 

% 

7' 

1. 19  x 1ft'  3 

4.03 

tion  of  the  aerodynamic  coefficient  time  history  estimates  reveals,  however,  that 
the  estimate  fit  is  better  in  the  longer  data  length  case.  C especially  exhibits 
the  best  fit  (Figure  6.30)  for  actual  time  history  data  that  lias  significantly  less 
well  behaved  variations.  The  coefficient  estimates  are  also  that  much  more  remark- 
able because  they  represent  estimates  of  two  dimensional  expansions  in  both  a and 
P with  simultaneous  variations  in  both  those  variables.  Two  affects  account  for  the 
improved  parameter  estimates: 

a.  A 50%  increase  in  data  length  represents  approximately  a similar 
increase  in  useful  information  for  identification  of  parameter  values. 

b.  Simultaneous  use  of  7 measurements  enables  use  of  what  information 
is  available  in  cross  mode  measurements.  This  effect,  although 
secondary,  does  contribute  to  the  ioial  information  something  not 
available  in  the  separated  lateral  and  longitudinal  mode  processing. 

6.5  PREDICTION 


The  ultimate  objective  of  all  parameter  identification  techniques  is  the 
ability  to  use  model  parameter  values  estimated  from  one  set  of  data  to  predict 
the  responses  of  a different  set  of  data.  Of  course,  the  operating  range  of  the 
data  to  be  predicted  should  not  extend  beyond  the  valid  limits  of  the  parameter 
estimates;  within  this  constraint,  the  predicted  response  should  closely  approxi- 
mate the  measured  response.  For  the  model  to  be  of  engineering  utility,  it  must 
not  only  be  able  to  reproduce  the  data  upon  which  it  is  based,  but  also  be  able 
to  closely  match  responses  to  different  inputs. 

Figures  6.31  through  6.34  show  the  predicted  and  actual  responses  of  the 
aircraft  using  parameter  estimates  from  lateral/longitudinal  response  case  1 for 
lateral/longitudinal  mode  case  2 (described  in  Figure  6.2).  The  inputs  are 


considerably  different  in  form  and  amplitude  between  cases  1 and  2 , which  can 
be  seen  by  comparing  the  time  histories  for  both  cases  . 
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Although  the  measurement  estimate  fits  for  the  prediction  are  not  quite 
as  good  as  the  original  parameter  estimation  data  fit,  the  correspondence  is, 
nevertheless,  more  than  satisfactory.  Even  the  aerodynamic  coefficient  time 
histories  reveal  remarkably  good  fits.  An  overall  assessment  of  this  relatively 
simple  31  parameter  estimate  model  must  be  that  it  is  successful  in  reproducing 
the  far  more  complex  simulation  of  real  aircraft  flight  characteristics. 

6 . 6 SUMMARY 

This  chapter  demonstrates  the  use  of  a new  procedure  for  the  identification 
of  stability  and  control  parameters  in  the  high  angle-of- attack  regime.  The  basis 
for  development  and  verification  of  this  procedure  has  been  noisy  data  generated  by 
a detailed  nonlinear  six  degree-of-freedom  simulator  for  the  F-4  aircraft.  This 
simulation  has  provided  a "controlled"  data  capability  which  allows  the  flexibility 
to  experiment  with  various  identification  problem  solutions . 

I 

The  method  has  been  applied  to  various  types  of  data  noise,  control  inputs, 
flight  conditions,  and  data  lengths.  Computation  has  been  directed  toward  expan- 
sion and  verification  of  the  capability  of  the  method,  and  detailed  quantitative 
assessments  of  various  trends  have  not,  within  the  scope  of  this  project,  been 
possible . 

a.  The  technique  for  isolating  and  estimating  values  of  significant  linear 
and  nonlinear  aircraft  model  parameters  from  realistic  flight  test  data 
has  been  developed  and  demonstrated.  The  subset  regression  step 
gives  the  overall  identification  procedure  the  unique  capability  of 
being  able  to  develop  parameter  estimates  without  any  a priori  know- 
ledge of  the  model  parameter  values  to  be  determined.  Although  some 
further  development  of  this  regression  technique  may  be  useful  to 
improve  structure  estimation  accuracy,  the  basic  principal  is  found 
to  significantly  improve  identification  estimates. 

b . The  procedure  is  robust  enough  to  yield  reliable  results  in  the 
presence  of  high  levels  of  noise,  both  measurement  and  process. 

c.  Only  significant  parameters  identifiable  from  the  given  data  are 
estimated.  Over-parameterization  is  thereby  avoided,  thus  producing 
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more  statistically  reliable  estimates,  as  well  as  eliminating  much 
extra  computation  that  would  produce  marginal  benefit  at  best. 

d.  Effects  of  input  design,  estimate  information  aggregation,  and  data 
length  are  shown  to  be  extremely  important  in  determining  the  informa- 
tion available  for  parameter  identification  in  a quantity  of  data.  For 
example,  input  design  results  indicate  that  small  amplitude  sinusoidal 

stabilator  deflections  give  good  C vs.  a non  linear  estimates.  These 

m 

estimates  are  further  improved  by  increased  data  length . 

e.  In  the  most  important  test  of  any  identification  technique,  the  ability 
of  the  procedure  to  develop  parameter  estimates  of  a quality  for 
excellent  model  predictions  is  demonstrated  . 
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VII.  CONCLUSIONS  AND  RECOMMENDATIONS 


b 


7 . 1 CONCLUSIONS 

7.1.1  Summary  of  the  Parameter  Identification  Procedure 


The  basic  objective  of  this  project  has  been  the  development  and  verification 
of  techniques  for  identification  of  the  stability  and  control  coefficients  of  high 
performance  aircraft  in  high  angle -of- attack  flight  regimes.  To  meet  this  objective 
the  following  three  tasks  have  been  completed: 

1 . Development  of  a detailed  digital  simulation  of  the  dynamic  response 
of  an  F-4  aircraft,  including  gust  and  instrumentation  noise  effects. 

2.  Application  of  an  advanced  regression  program  to  determine  the 
significant  aerodynamic  contributions  to  high  angle- of- attack  re- 
sponses . 


Implementation  of  a highly  versatile  version  of  the  maximum  likeli- 
hood parameter  identification  algorithm  which  estimates  the  values 
of  significant  stability  and  control  coefficients  for  the  high  angle-of- 
attack  flight  regime. 


1 he  accomplishment  of  these  tasks  ha  demonstrated  that  significant  aero- 
dynamic contributions  to  complicated  aircraft  responses  can  be  estimated  by  the 
sequential  process  of:  (a)  determining  which  effects  are  important  (e.g.  , the 
i egression  phase),  and  (b)  determining  highly  accurate  parameter  estimates  of 
these  effects  (e.g  . , the  identification  phase)  . The  importance  of  this  new  combined 
approach  is  that  the  source  of  many  obstacles  to  identifying  high  angle-of- attack 
coefficients  from  flight  data  is  avoided— namely , problems  of  identifiabilily  and 
over-parameterization.  Ibis  approach  is  no.  only  useful  for  producing  rapid 
determination  of  significant  causes  of  complicated  urcraft  response , but  also  for 
identifying  the  parameters  of  a wider  class  of  nonlinear  systems  than  has  previous- 
ly been  reported . 


I 
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7.1.2  Specific  Chai-acteristics  of  the  Integrated  Process 


Application  of  the  model  estimation  program  demonstrates  a significant 
improvement  in  maximum  likelihood  efficiency,  both  in  required  computation  and 
in  accuracy  of  the  results . The  improvement  in  performance  of  the  entire  identi- 
fication process  is  based  on  the  following  characteristics  of  the  regression 
method: 


The  subset  regression  method  selects  parameters  on  their  ability  to 
match  the  measured  response.  By  selecting  the  minimum  number  of 
variables  to  accomplish  this  match,  the  most  significant  variables 
are  isolated  . 


The  method  yields  a px-iori  estimates  from  the  data  above,  and  does 
not  itself  require  initial  estimates.  Though  usually  biased,  the  esti 
mates  given  by  the  program  are  frequently  better  than  a priori 
estimates  from  other  sources. 


The  method  gives  significant  evaluation  of  tire  estimated  parameters, 
which  serves  an  essential  role  in  the  final  selection  of  a model  struc 
ture  to  be  used  in  tire  maximum  likelihood  algorithm. 


The  computational  requirements  of  the  regression  algorithm  are 
modest,  and  the  algorithm  can  be  used  quickly  to  not  only  evaluate 
input  designs,  but  also  effects  of  measurement  errors  . 


Subsequent  applications  of  the  maximum  likelihood  algorithm  are  optimized 
by  introduction  of  the  regression  "pre-processing"  . The  implementation  of  the 
maximum  likelihood  for  this  work  can  accept  either  linear  or  nonlinear  aerody- 
namic expansions  of  the.  force  and  moment  coefficients.  Further,  this  algorithm 
can  process  parameters  for  only  one  dynamic  equation  or  all  six.  The  program 
includes  calculations  of  parameter  significance  level. 
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7 • 1 • 3 Application  to  Simulated  Data  for  the  High  Angle-of-Attack  Stall/Post- 
Stall  Regime 

The  evaluation  of  the  integrated  parameter  identification  process  has  been 
conducted  with  simulated  data.  The  simulation  is  highly  nonlinear  and  of  much 
higher  order  than  the  identification  models  used.  The  measurement  model  is  also 
highly  nonlinear  and  not  readily  describable  by  linear  error  models  for  the.  entire 
flight  regimes  tested.  Behavior  of  its  modeling  detail  output  of  the  simulation  is 
regarded  as  equivalent  to  that  obtained  from  actual  flight  test  data,  at  least  for  the 
purposes  of  the  identification  evaluation. 


The  main  conclusion  of  applying  the  identification  procedure  to  the  data  is 
the  predominant  requirement  for  linear  aerodynamic  models,  even  for  "nonlinear" 
regimes.  Recognizing  that  such  conditions  may  not  be  attainable,  nonlinear 
models  are  required  for  large  inputs  where  angle-of- attack  and  sideslip  are  simul- 
taneously excited,  such  as  would  be  obtained  for  large  amplitude  sinusoidal  inputs 
or  special  maneuvers . The  regression  step  is  essential  since  it  indicates  whether 
or  not  nonlinear  contributions  are  required  . The  ability  to  recognize  and  identify 
such  nonlinear  terms  is  requisite  for  determination  of  coefficients  such  as  C 

m 

and  C . 
n 

Primary  emphasis  has  been  placed  on  the  identification  of  static  force  and 
moment  coefficients.  Greatest  errors  have  been  found  in  estimates  of  the  dynamic 
rate  derivatives  and  control  effectiveness  derivatives . These  errors  are  attributed 
to  the  relatively  low  rates  which  were  obtained  for  inputs  to  identify  the  important 
static  derivatives.  Inputs  may  be  modified,  using  input  design  algorithms  dis- 
cussed in  the  report,  to  improve  the  estimation  of  these  derivatives.  It  must  be 
noted  that  the  errors  of  the  dynamic  derivatives  which  were  obtained  could  not 
be  found  by  comparing  time  history  matches  of  the  measurement  of  rate,  as  these 
were  excellent  in  all  cases.  Only  with  comparison  of  the  actual  parameter  time 
histories  were  these  errors  observed. 

Techniques  for  improving  parameter  estimates  include  input  design,  in- 
creased data  length  (where  possible),  and  data  aggregation  (e.g.,  statistically 
combining  results  of  several  tests)  . Such  modifications  have  a significant  effect 
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on  flight  test  planning  and  should  be  available  for  identification  of  important  but 
difficult  to  excite  parameters.  The  integrated  parameter  identification  pi'ocess 
as  developed  in  this  work  is  amenable  to  such  flight  test  design,  using  procedures 
detailed  in  the  previous  chapters. 

The  identification  procedure  is  robust  with  respect  to  process  and  measure- 
ment noise.  For  this  reason,  the  option  of  identifying  process  noise  statistics  was 
not  included  in  the  maximum  likelihood  implementation  (and,  hence,  a substantial 
reduction  in  computer  time)  . 

The  confidence  established  in  the  procedure  led  to  evaluation  of  the  identi- 
fication results  by  a prediction  criterion.  Specifically,  the  parameter  estimates 
from  a specific  lateral  maneuver  with  one  input  are  used  to  predict  the  response 
of  the  simulation  to  another  input  (for  roughly  the  same  flight  regime)  . This 
prediction  capability  was  verified  with  excellent  results , 

In  summary,  it  is  concluded  that  the  application  of  the  integrated  parameter 
identification  process  developed  for  the  high  angle- of- attack  stall/post-stall 
regime  offers  significant  improvements  in  the  ability  to  identify  not  only  para- 
meters, but  also  the  entire  system  structure  and  parameters. 

7 . 2 RECOMMENDATIONS 

The  method  developed  under  this  contract  promises  to  be  a valuable  tool 
for  Navy  flight  test  data  reduction  requirements.  As  discussed,  the  procedure  is 
highly  oriented  toward  user  needs  and  applicable  to  both  linear  and  nonlinear 
models  of  typical  Navy  aircraft.  This  versatility  is  essential  for  many  operations, 
tests,  and  evaluation  requirements  for  high  performance  aircraft.  Typical  applica- 
tions, for  example,  include: 

1.  Documentation  of  aircraft  stability  and  control  characteristics  for 
critical  flight  regimes  (e.g.  , high  angle- of- attack)  or  operations 
(e.g.  , carrier  landing).  Such  documentation  is  required  fjr  design 
evaluation,  transfer  function  formulation  (useful  in  handling,  quality 
evaluation),  and  maneuver  envelope  expansion. 
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2.  Determination  of  effects  of  special  design  modifications  such  as  lead- 
ing edge  slats,  advanced  avionics  control  systems,  or  modified  air- 
craft instrumentation.  The  method  developed  for  this  project  will  be 
applicable  for  investigating  and  evaluating  such  modifications  which 
are  incorporated  to  improve  maneuverability  or  handling  qualities. 

The  method  and  analyst's  discussed  in  this  report  is  also  capable  of  being 
expanded  to  particular  objectives.  Such  objectives  include: 

1.  Combination  of  the  regression  algorithm  with  the  maximum  likelihood 
program  to  give  one  single  program.  The  requirement  here  is  to 
design  the  interface  between  the  two  techniques  such  that  the 
amount  of  "tuning"  the  two  methods  to  each  other  is  included  auto- 
matically . 

2.  Approximation  analysis  to  the  maximum  likelihood  algorithm  to  allow 
on-line  and/or  real  time  capability  to  provide  faster  acquisition  of 
parameters . The  requirement  is  to  more  completely  model  the  data 
links  of  on-board  or  ground-based  measurement  systems  and  to  modify 
the  maximum  likelihood  algorithm  to  require  less  computer  capacity 
while  operating  in  real  time. 

3.  Continuation  of  input  design  development,  for  nonlinear  systems  to 
reduce  flight  test  time  as  well  as  improve  parameter  estimates.  The 
requirement  is  to  relate  the  orthogonal  polynomial  representation 

to  meaningful  aerodynamic  flight  test  terms. 


4 . Implementation  of  maximum  likelihood  techniques  for  fault  detection 
of  multivariable  systems.  The  requirement  is  classification  of  faults 
and  their  relationship  to  parameters  of  the  system. 

5.  Optimization  of  control  surface  and  instrument  location.  The  maxi- 
mum likelihood  technique  uses  several  functions  which  depend  on 
control  surface  locations  and  instrument  location.  Optimization  of 
these:  functions  will  lead  to  strategies  for  not  only  improving  the  iden- 
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tifi cation  process,  but  also  for  the  general  control  synthesis  problem 
of  light  weight  fighters,  for  example.  Such  aircraft  must  make  opti- 
mum use  of  control  surfaces  for  satisfactory  maneuverability. 

Modification  of  the  models  for  aircraft  could  easily  be  made  to  extend 
the  integrated  parameter  identification  to  missiles,  RPV , or  VTOL 
aircraft . 

Other  vehicles  which  can  be  treated,  without  extensive  development, 
are  surface  ships  (conventional,  surface  effect  ships,  semi-submerged 
hull  designs,  and  hydrofoils)  as  well  as  submarines.  Such  systems 
are  well  known  to  be  susceptible  to  significant  nonlinear  hydrody- 
namic nonlinear  forces  and  moments. 

One  promising  application  of  identification  is  the  determination  of 
aircraft  engine  parameters.  The  high  performance  requirements  of 
such  engines  is  placing  stringent  accuracy  requirements  on  their 
design  and  test.  The  potentials  of  engine  multi- variable  control 
augment  these  accuracy  goals.  Modeling  and  identification  of  turbine 
dynamics,  aerodynamics,  and  thermodynamics  could  greatly  aid  in 
achieving  these  requirements. 

I 


: 

i 
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APPENDIX  A 

MAXIMUM  LIKELIHOOD  IDENTIFICATION  OF  PARAMETERS 


A . 1 INTRODUCTION 

The  maximum  likelihood  method  is  one  of  the  most  flexible  techniques 
in  statistics  for  identification  of  parameters  from  input-output  data.  Suppose 
it  is  possible  to  make  a set  of  observations  on  a system,  whose  model  has 
p unknown  parameters  0.  For  any  given  set  of  values  of  the  parameters  G 
from  the  feasible  set  0,  we  can  assign  a probability  p(Z|8)  to  each  outcome 
Z.  If  the  outcome  of  an  actual  experiment  is  z,  it  is  of  interest  to  know  which 
sets  of  values  of  0 might  have  led  to  these  observations.  This  concept  is 
embedded  in  the  likelihood  function  .'/(0|z),  This  function  is  of  fundamental 
importance  in  estimation  theory  because  of  the  likelihood  principle  of  Fisher 
and  others  [40-42]  which  states  that  if  the  system  model  is  correct,  all  informa- 
tion about  unknown  parameters  is  contained  in  the  likelihood  function.  The 

A 

maximum  likelihood  method  finds  a set  of  parameters  0 to  maximize  this  likeli- 
hood function 


0 = max '/X  0/z) 
OcO 


(A.l) 


In  other  words,  the  probability  of  the  outcome  z is  higher  with  parameters 

A 

0 in  the  model  than  with  any  other  values  of  parameters  from  the  feasible 
set.  Usually  it  is  more  convenient  to  work  with  the  logarithm  of  the  likelihood 
function  (it  is  possible  to  do  so  because  the  logarithm  is  a strictly  monotonic 
function)  . 


The  great  asset  of  the  maximum  likelihood  method  is  that  it  can  be  used 
with  linear  or  nonlinear  models  in  the  presence  of  process  and  measurement 
noise.  The  maximum  likelihood  estimates  arc  asymptotically  unbiased, 
consistent  and  efficient. 

A.  2 GENERAL  NONLINEAR  SYSTEMS 


Consider  the  general  nonlinear  aircraft  equations  of  motion 


^v""  !S*v  >*r«v  - ,'• 


. ..  


x = f(x  , u , 0 , t)  + r(0,t)s 


0 < t < T 


E(x(0))  = xQ(0) 


E{(x(0)  - x (0))(x(0)  - x (0))r)  = P (0) 
° o o 


where 


x(t)  is  n x 1 state  vector 


u(t)  is  5 x 1 input  vector 


9 is  p x 1 vector  of  unknown  parameters 


a is  n x q process  noise  distribution  matrix 


w(t)  is  q x 1 random  process  noise  vector 


Sets  of  m measurements  y(t^)  are  taken  at  discrete  times  t^_ 


y(t,r)  = h(x(tk),u(tk),  0,tk)  + v(t,.) 


k = 1,2,3 N 


(A.  3) 


w(t)  and  v(t  ) are  Gaussian  random  noises  with  the  following  properties 


E(w(t))  = 0 E(v(tk))  = 0 E(w(t)  vT(t  ))  = 0 


rp  r 

E(w(t)  w (x))  = 0(0, t)  5(t-x)  , E(v  (t.)  vT(t,  ))  = R (0  ,t.  )5  , (A. 4) 

J k J jk 


The  unknown  parameters  are  supposed  to  occur  in  the  functions  f and 


h and  in  matrices  T,  Q,  R,  Po  and  xq.  In  the  analysis  to  follow,  the  model 


and  the  functional  form  of  f and  h is  assumed  known  correctly. 


The  set  of  observations  yCt^.yCtp yCt^)  constitutes  the  outcome 

Z m this  case.  The  likelihood  function  for  0,  which  has  the  same  form  as  the 
probability  of  the  outcome  z for  a certain  value  of  parameters  0,  is  given  by 


s?(e|z)  s p(z|0) 

= pCyCtj).  y(t2), . . . ,y(tN)|o) 

= p(yn|0)  (A. 5) 

where 


Yk  = k = i > 2, . . . ,n 

p(yn|0)  = p(y(tN)|YN.1,0)-P(YN_1|0) 

p(y(tN)|YN-i*0)  P(y(tN-l)lYN-2'0) 

N 

= 17  p(y(t.)|Y  .,0)  (A. 6) 

i=l  1 1 

The  log-likelihood  function  is 


N 

log  (5?(0|z))  = Z l0g'{p(y(ti)|Y._1.e)}  + constant  (A. 7) 

i=l 

To  find  the  probability  distribution  of  y(t.)  given  Y^_ ^ and  0,  the  mean  value 
and  covariance  are  determined  first. 

E(y(t.)|Y._1(e)  ^ y (i/i-1)  (A. 8) 

The  expected  value  or  the  mean  is  the  best  possible  estimate  of  measurements  at 
a point  given  the  measurements  until  the  previous  point. 


cov 


(y(ti)lYi_l>e)=  E{  (y(t.)-y (i/i-l))  (y(t,)-y (i/i-l))T) 


A E { V (i ) v(i)T} 


(A.  9) 
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A B(i) 


v(i)  are  the  innovations  at  point  i and  B(i)  is  the  innovations  covariance. 
Since 


y(t.)  - E(y(ti)|Y._1,0)  = v(i) 


(A. 30) 


it  follows  that  y(t.)  given  Y. and  0 have  the  same  distribution  as  V (i) . 
Kailath  [44  ] has  shown  uiat  as  the  sampling  rate  is  increased,  the  innovations 
v(i)  tend  towards  having  a Gaussian  density.  Assuming  a sufficiently  high 
sampling  rate,  the  distribution  of  v(i)  and,  therefore,  y(h)  given  Y. and  G 
is  Gaussian , i . e . , 


pCy(ti)|Y._1,e)  = 


exp|-  Iv(i)TB  ](i)  v(i)j 
(27t)m/2|B(i)|J/2 


(A.  11) 


log  {P(yCti)|Yi  l,G)}  = - I v(i)TB_1(i)v(i)-i  logjB (i)|  + constant 

(A.  12) 


The  log-likelihood  function  of  Equation  (A.  7)  can  be  written  as 


log  (^(6|z)) 


\ 2 {vT(i)  B'](i)  v(i)  + logjB (i)|) 
i=l 


(A.  13) 


An  estimate  of  the  unknown  parameters  is  obtained  by  maximizing  the  like- 
lihood function  or  the  log-likelihood  function  in  r feasible  set  of  para- 

meter values . 


G = max  log( '/'(G|z)) 
0c0 


(A. 14) 


= max  - \ E {vT(i)B  1 (i) v (i)  + logjB(i))} 
Ge©  l 6 i=l 


(A.  15) 


The  log-likelihood  function  depends  on  the  innovations  and  their  covar- 
iance. To  optimize  the  likelihood  function,  a way  must  be  found  for  determin 
ing  these  quantities.  Both  innovations  and  their  covariance  are  outputs  of  an 
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The  extended  Kalman  filter  is  conventionally  divided  into  two  parts  . In 
the  first  part,  called  the  prediction  equations,  the  state  equations  and  state 
estimate  covariance  equations  are  propagated  in  time  from  one  measurement 
point  to  the  next.  In  the  second  part,  called  the  measurement  update  equations, 
the  measurements  and  associated  measurement  noise  covariances  are  used  to 
improve  state  estimates.  The  covariance  matrix  is  also  updated  at  this  point 
to  reflect  the  additional  information  obtained  from  the  measurements. 

Prediction  Equations 

The  state  prediction  is  done  using  the  equations  of  motion  (A. 2). 

Starting  at  time  h_j  with  current  estimate  x(i- lji- 1)  of  tire  slate  x(t.)  anu  the 
covariance  P(i-l|i-l),  the  following  equations  are  used  to  find  the  predicted 
state  x (i| i-1)  and  the  associated  covariance  P(i|i-T);  see  Bryson  and  Ho  [45]. 


dt  = fC>Kt|  ti_1)  ,u(t) , 0,t> 


(A.  16) 


P^ti-i)  = F(t)  P(t|t._1)  + PCt)t._1)  FT(t)  + ror"- 


(A.  17) 


t.  , < t < t. 
i-I  - - i 


The  nxn  matrix  F is  obtained  by  linearizing  f about  the  best  current 


estimate 


F (t)  = 


9f(^(t|tj_]),u(t),0,t) 


3x(t|t._1) 


(A.  18) 


Using  (A.  16)  to  (A.  18),  we  can  obtain 


xU.lt._j)  A x(i| i-1) 


and  PU.lt._j)  A P(i| i-1) 


(A. 19) 
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Thereafter,  the  measurement  update  equations  are  used 


Measurement  Update  Equation 


The  covariance  and  state  estimate  are  updated  using  the  measurements. 
The  necessary  relations  are  derived  by  Bryson  and  Ho  [45  ] and  are  presented 
here  without  proof.  The  innovation  and  its  covariance  are 


where  H is  obtained  by  linearizing  h 


The  Kalman  gain  and  the  state  update  equations  are 


Finally,  P(i|i),  the  covariance  of  error  in  updated  state,  is  obtained  by 


One  can  now  return  to  the  time  update  or  prediction  equations 


A.  3 OPTIMIZATION  PROCEDURE 


Many  possible  numerical  procedures  can  be  used  for  this  optimization  prob 
lem . Modified  Newton -Raphson  [20]  or  Quasilinearization  [21]  have  been  found 
by  experience  to  give  quicker  convergence  than  most  procedures  like  the  con- 
jugate gradient  or  the  Davison  method.  The  modified  Newton-Raphson  is  a 


second  order  gradient  procedure  requiring  computation  of  first  and  second 
order  partials  of  the  log-likelihood  function. 


The  last  three  terms  in  the  equation  for  second  partial  of  the  log-likelihood 
function  involve  second  partials  of  innovation  and  its  covariance.  Those  terms 
are  usually  dropped.  So  the  second  partial  is  approximated  by 


a log  (tf(eiz))  - 


ae.  ao, 

i k 


^ avr(i)n-i,.N  av(i)  avT(i)  „-i,.x  dECO^-i,^  .,.N 

,Mir-B  (l)  ”a-o---aG7-B  (l)^o-B  (l)v(l) 

i=1  j k k ] 


a-m  B-ki)  v(o 

l k 


+ vTCi)B_1(i)  B_1(i)  ^—-B  JCi)  v(i) 

i k 


1 m „-l,.„  3B(i)  „-l,..x  8B(i) 

' 2 Tr  B (l)  ~ae~  B (0  “aor  ( 
L J kj 


(A.  28) 


The  gradients  of  innovation  and  its  covariance  for  parameter  are  (Equations 


(A. 20)  and  (A. 21)) 


3v(i)  3h  | 


8x(i|i-l)  3h 


0XI  x=x(i[i-l) 


(A. 29) 


j = 1,2, ... ,p  ; i = 1,2, . 


= p(i|i-i)HTa)  + mo  0B8(l^-1)  HT(i) 
j j j 


. TT/.x  T,,.|.  ^ ai-iT(i)  , 3R 
+ H(i)  P(i|i-1)  -QQ—  + §Q- 


(A.  30) 


j = 1,2, ...  ,p 


i = 1,2,  . . . ,N 


Recursive  equations  can  be  obtained  for  gradients  of  the  predicted  state 
and  its  covariance.  This  is  done  in  stages  by  using  the  prediction  and 
measurement  update  equations  of  Kalman  filter.  Differentiating  (A.  16)  to  (A.  18) 


with  respect  to  0. 
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! 


d 8x(t|t.j)  _ af(x(t|t._1),u(t),t,o)  af(x(tlt._]),u(t),t,e) 

~ + 


dt  ao. 

j 


axait.^) 


9x(t|t._1) 


ao. 

3 


9x(0|0) 

ao. 

3 


dx  ( 0 ) 
o 

ae. 

3 


(A. 33) 


d ap('IVi> 

dt  30. 


OFW  rrtl  , , p.n  rT . , 

-go—  + F(l)  uW. + — aa F «> 

3 3 3 


+ i1  QrT  + + ^ 


3 3 


dr1 

ao. 


0P(O|O)  _ aPo(0) 

ao.  ae. 

. 3 3 


t.  , < t < t 
1-1  - — 


(A. 32) 


9F(t)  _ aZf(x(t|t._1),u(t),  9,0 


ao. 

3 


ao.  dx  (tjt.  , ) 


3 


-V 


j = 1,2, ...  ,p 


(A. 33) 


The  sensitivity  functions  are  updated  at  measurement  points  by  differentiating 
(A.  22)  to  (A.  25)  with  respect  to  0^ . 


9H(i) 

ao. 

3 


a2h(x(i|i-l),u(t.),0,t.) 
ao.  ax(i[i-i) 


(A.  34) 


3K(i) 

ao. 

3 


DP(i|i-l) 

DO. 

3 


HT(i)  B~ 1 (i ) + P(i|i-1)  B-1(i) 

3 


-PCifi-])  H T (i ) rT-’o)  B_1(i) 

j 


(A.  35) 


177 


The  negative  of  the  matrix  of  second  partials  of  the  log- likelihood  func 
tion  is  called  the  information  matrix  M.  The  step  size  A0  for  parameter  esti- 
mates is  given  by 


The  information  matrix  provides  a lower  bound  on  parameter  estimate  covar 
iances  , i .e . , 


This  is  called  the  Cramer-Rao  lower  bound . 
estimates  approach  this  bound  asymptotically 


The  maximum  likelihood 


A.  4 LINEAR  SYSTEMS 


In  a li  near  system , the  functions  f and  h are  defined  as 


In  this  case,  it  has  been  shown  that  if  the  model  is  correct,  the  innovations 
are  white  and  have  Gaussian  density  at  the  true  values  of  the  parameters. 
Therefore,  the  assumption  of  fast  sampling  rate  is  not  necessary. 


The  basic  algorithm  is  the  same.  However,  some  of  the  equations  can 
now  be  simplified.  The  equivalence  between  F(t)  and  H(i)  of  (A.  18)  and 
(A.  22)  and  F(0,t)  and  H(0,t.)  of  (A.  40)  is  obvious.  Equations  (A.  16)  and 
(A. 20)  now  become 


Equations  (A. 29)  and  (A. 31)  can  be  written  as 


All  other  equations  remain  the  same 


There  is  considerable  reduction  in  computation  requirement  for  time- 
invariant  linear  system.  In  this  case,  matrices  F,  G,  II,  D,  F , Q and  R and 


their  derivatives  with  respect  to  parameters  are  constant.  Gupta  [ 48  ] has  shown 
that  the  computation  of  state  sensitivities  can  be  reduced  to  many  fewer  equa- 
tions. Similar  reductions  are  possible  in  the  computation  of  covariance 
sensitivities . 


A.  5 TIME  INVARIANT  LINEAR 


Y STEMS  IN  STATISTICAL  STEADY  STATE 


In  many  aircraft  applications , the  Kalman  filter  is  in  steady  state  for  the 
duration  of  the  experiment.  This  occurs  when  the  Kalman  filter  is  in  opera- 
tion for  a sufficiently  long  time  and  the  process  and  measurement  noise  covari- 
ances do  not  change.  The  Kalman  gain  and  the  innovations  and  the  state 
covariances  approach  constant  values.  The  time  update  and  measurement  up- 
date equations  for  the  covariances  are 


By  definition  of  the  steady  state 


Therefore 


I ' 


A cp(At)  (I-KJj)  P(i|i-1)  cp( At)  + Q' 


(A.  49) 


= cp (At)  (P(i|i-1)  - KBK  ) 9 (At)  + Q' 


Using  (A.  46),  (A.  47)  and  (A.  49),  we  can  solve  for  P(i|i-1)  and  then  find 
K and  B.  Also,  it  can  be  shown  that 


3P 


= A. 


9P  , T 


A,  + A_  - cpPA.Pcp ' 


ae.  i ae.  “i  ‘*2  ^ “3" 

3 3 


CA.50) 


= (l-KH)  HT(HPHT  + R)  1 + A2 


3 


(A. 51) 


3 


= phT  + H dP_  „T  3HT  3R 

ae.  ae.  + H air  11  + Hr  ao7  + ae: 

3 j j 3 j 


(A.  52) 


where 


Ax  = <p(I-KH) 


T 


A2  " le~  (I"KH)Pcf)I  + <P(I-KH)P  §J-  - (pK  ~ PcpT  + |Q 


3 


3 


ae. 

3 


a3  = He  b_1h  + hTb_1  sir  P,,T  + HP  Hr  + STT  B_]h  (A. 53) 


3 


ae. 

3 


ae.  ae. 

3 3 


P A P(iji-l) 


Thus,  it  is  possible  to  solve  for  using  (A.  50)  and  then  find  4^-  and  ~ 

ou-  00.  00 
J »\p  j j 

from  (A. 51)  and  (A. 52),  (A. 50)  is  a linear  equation  in  t-tt-  and  the 

dO. 

1 


coefficient  of  the  unknown  matrix  does  not  depend  on  the  parameter  6 . Thus 

3 

the  sensitivity  of  stale  covariance  matrix  can  be  determined  very  quickly  for 
all  parameters . Once  the  sensitivity  of  P , K and  13  for  unknown  parameters 
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is  determined,  only  state  sensitivity  equations  need  to  be  updated.  The  compu- 
tation of  state  sensitivity  functions  can  be  reduced  to  many  fewer  equations  as 
pointed  out  in  Section  A. 4. 

An  approximation  suggested  by  Mehra  [ 34  ] simplifies  the  problem  further 
The  unknown  parameters  are  defined  to  include  elements  in  K and  B matrices 
instead  of  Q and  R.  Optimizing  the  log-likelihood  function  for  parameters  in 
B gives 


N 

I v(i) 
i=l 


v1  (i) 


(A . 54) 


Ihe  gradient  of  the  log-likelihood  function  with  respect  to  other  unknown 
parameters  is 


9 Log  ( '/'(6[  ?,)) 

ae. 

i 


’ ry>  A 

-L  v (i)  B 
i=l 


9v(i) 

ae. 

3 


(A. 55) 


The  sensitivity  of  innovations  to  parameters  is  determined  using  the  following 
recursive  equations 


dt 


x(t|t._1)  = Fx(t|ti_1)  + Gu(t) 


d f)F  9x(t|t.  ) .... 

d t ~ bit:  90:  + F ~9o: — + m~  u(t) 

3 J J j 


j = 1,2,  . . . ,p  t.  . < t < t. 

l-l  — — i 


(A.  56) 


v(i)  = y(i)  - Hx(i|i-1)  - Du(t.) 
x(i|i)  = x(i|i-l)  + Kv(j) 


(A. 57) 
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Sv  (i)  bn_ 

ae.  " ao. 

3 J 


x(i|i-l) 


be. 

j 


9D 

80. 

J 


u(t.) 


3*010  _ ax Ci|i~  1 ) bK_  rn  , K bv(i) 

be.  bn.  be.  be. 

3 J J 3 


(A. 58) 


Note  that 


9K  0 . 

KTi-  =0  if  0.  is  net 

bO.  j 


element  cf  K matrix 


if  0.  A K. 

3 “ 3 k 


(A. 59) 


where  I. . is  a matrix  of  all  zeroes  except  a 1 at  the 
3 k y 


jOk'  position. 


This  approximation  simplifies  the  optimization  considerably.  However, 
this  usually  leads  to  an  overparameterized  model.  In  other  words,  once  K and 
3 are  determined,  it  is  not  possible  to  find  any  corresponding  T,  Q and  R 
which  have  the  desired  (known  a priori)  structure.  A good  estimate  of  ele- 
ments in  P,  Q and  R matrices  can  be  obtained  by  a least-squares  type  approach. 
The  fit  to  the  obsei*ved  data  is  better  than  with  true  values  of  parameters 
but  the  parameter  estimates  do  not  have  minimum  variance.  This  approximation 
is  not  good  in  aircraft  application  where  the  structure  of  Q and  R is  known 
fairly  well,  but  is  excellent  where  there  are  many  process  noise  sources  and 
the  characteristics  of  both  Q and  R are  relatively  unknown  (e.g.,  economic 
systems)  . 


A. 6 MAXIMUM  LIKELIHOOD  WITH  NO  PROCESS  OR  MEASUREMENT  NOISE 


Tile  maximum  likelihood  method  can  be  simplified  when  cither  process 
noise  or  measurement  noise  are  absent. 


...... 


No  Process  Noise 


If  the  process  noise  is  zero  and  initial  states  are  known  perfectly,  i.e., 
w(t)  and  P(0)  are  zero,  the  covariance  of  the  error  in  the  predicted  state 
is  also  zero.  It  is  clear  from  (A. 23)  that  Kalman  gains  are  zero.  The  inno- 
vations is  the  output  error,  i.e., 


v (i)  = y 0)  - h(x(t.),u(t.),G,t.) 


and  the  innovation  covariance  is  (A.  21) 


B(i)  = R 


the  log-likelihood  function  is , 


(A.  60) 


(A. 61) 


N 


log  mo |z))  = 4 S vT(i)  R_1v(i)  + log  |R| 
* i=l 


which  on  optimizing  for  unknown  parameters  in  R gives 


N 


(A. 62) 


R=  h r v(i)  vT(i) 

N i=l 


(A. 63) 


The  equality  in  (A. 63)  holds  only  for  those  elements  of  R which  are  not 
known  a priori.  For  instance,  even  if  R is  known  to  be  diagonal,  the  right 
hand  side  matrix  will  not  be  diagonal  in  general,  but  the  off-diagonal  terms 
should  be  ignored  before  they  are  equated  to  R.  Using  (A. 63)  in  (A. 62) 


1 N T -_i 

log  ( ,^(0|z))  - ” y T v (i)  R v(i)  + constant 

6 i=l 


(A. 64) 


The  optimizing  function  is  the  same  as  that  for  the  output  error  method  except 
that  the  measurement  noise  covariance  matrix  is  determined  using  (A. 63) 
and  is  used  as  the  weighting  matrix  in  the  criterion  function.  In  the  output 
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error  method,  the  measurement  noise  is  assumed  known  and  the  weighting 
function  is  arbitrary. 

The  first  and  second  derivatives  of  the  log-likelihood  function  with 
respect  to  unknown  parameters  are 


g§-log(tf(0|z))  = - S vT(i)  R"1 


(A.65) 


9 log(.ff(0[z)) 

90.  90, 

J k 


v j9vT(i)  A-1  9v(i) 


* *T<» n SH- 

J k 


The  terms  in  the  second  derivative  are  approximated  as 


9 log(<y(0[z)) 

90.  90, 

J k 


- E |9vTW  r"1  9v(i)l 

i=l(  d9k  80j  ) 


(A.  66) 


(A. 67) 


No  Measurement  Noise 

If  all  states  are  measured  with  no  noise,  the  covariance  of  the  error  in 
state  estimates  is  zero  at  the  beginning  of  any  time  update. 


P(i-l|i-l)  = 0 


and  x(i-l|i-l)  = x(i-l) 


(A. 68) 


It  is  easy  to  show  in  this  case  that  for  fast  sampling  the  log-likelihood  function 
is  quadratic  in  the  difference  between  measured  values  of  x and  f(x  ,u,0,t) . 

The  method  reduces  to  the  equation  error  method , the  weight  W being  chosen 


• Vi  . 
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Thus,  the  maximum  likelihood  method  and  equation- error  methods  are  equiva 
lent  except  for  the  technique  for  choosing  the  weighting  matrix. 


APPENDIX  B 

MODEL  STRUCTURE  DETERMINATION  BY  MULTIPLE  REGRESSION  TECHNIQUES 


This  appendix  provides  an  intuitive  and  mathematical  basis  for  the 
multiple  regression  technique  application  to  model  building  from  data.  The 
basic  reference  is  Volume  II  of  Kendall  and  Stuart  [49]. 


B.l  LEAST  SQUARES  METHOD 

The  method  of  least  squares  (also  known  as  regression  method  or  equa- 
tion error  method)  is  based  on  the  minimization  of  the  scalar  sum  of  squares 


J = (y  - X0)  (y  - X0) 


(B.l) 


with  respect  to  the  components  of  0.  The  cost  (B.l)  is  predicated  on  a 
linear  model  of  the  form 


y = X0  + e 


(B  .2) 


where  y is  an  (mxl)  vector  of  observations,  X is  an  mxp  matrix  of  known 

coefficients*  (m  > p)  and  0 is  a pxl  vector  of  parameters,  e is  an  (nxl) 

vector  of  disturbance  or  error  random  variables  with  <o(e)  = 0 (zero  mean) 
T 2 

and  variance  <?(ee  ) = o I (components  of  e are  uncorrelated  with  the  same 
variance  o^)  . 


Minimization  of  J with  respect  to  0 yields  the  well-known  least  squares 
estimator 


a rF  —IT* 

0 = (XXX)  X y 


(B . 3) 


The  discussions  of  this  appendix  assume  the  data  and  parameters  are 
centered  about  their  means . 
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I 


■ rp 


where  (X  X)  is  nonsir^gular . Substituting  Eq.  (B.2)  into  Eq.  (B.3),  it  is 
easily  shown  that 


<?(0)  = e 


<?{(0-0)(0-e)T}  = a2(XTX)"1 


(B . 4) 
(B . 5) 


where  X is  fixed  in  repeated  samples  and  assuming  <?(X  e)  = 0. 


2 . 


In  actual  experimentation,  a is  not  known  and  must  be  estimated. 

This  is  due  to  the  fact  that  e = y - X0,  the  error,  is  a nonobservable,  stochas- 


tic variable.  However,  estimates  of  the  variance  of  e = y - X0  (Fig.  B.l).  s , 


may  be  determined  which  estimate  a . The  sum  of  squares  of  the  residuals 


aT  a 

is  found  by  solving  Eqs.  (B.l)  and  (B.2)  for  e e,  i.e., 


eTe  = (y  - X0)T(y  - X0)  = 


(B . 6) 


T 

= y 


where  ,Ji-  (I  - X(X^X)  *X^)  . Note  that  --yff  - . Since  <A(  is  thus 


m 


idempotent  and  \Xi4(  - m - p,*  taking  the  expectation  of  Eq.  (B.6)  yields 


21  A T * *Ta 

s*  =~±—  (y-X0)1(y-X0)  = ^ 

m-p  3 3 m-p 


(B . 7) 


which  estimates  o . Equation  (B.7)  in  Eq.  (B.5)  gives  the  sample  estimate 

A 

of  the  covariance  of  0. 


(0-0)  (0-0)T  = s2(XTX)_1 


(B  .8) 


1 


A 

Note  that  Eq.  (B.3)  produces  unbiased  estimates  of  0 only  if  the  model  is 
correct . 

* T -IT 

The  rank  of  Iv/K=  X(X  X)  X is  p,  and  the  rank  of  Im  is  m . If  the 
observations  are  not  centered,  then  p is  replaced  by  p+l- 
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B.2  DECOMPOSITION  OF  VARIANCE 


Substituting  the  estimate  0 into  J,  it  is  found  that 


y aT  A T aT  A 

J = e e = y y - y y 


(B  .9) 


where 


y = (y|X)  = X6  (B.10) 

Equation  (B.9),  a consequence  of  the  orthogonality  theorem  of  least  squares, 

T 

decomposes  the  sum  of  squares  of  the  data  variation  (y  y)  into  the  contribu- 

aT  a 

tion  from  the  sum  of  squares  of  the  regression  equation  (y  y)  and  sum  of 
squares  of  the  residuals  (e^e),  i.e., 

yy  = yy+ee  (B.  11) 

The  basic  idea  of  subset  regression  is  to  compute  the  reduction  in 

aT  a 

error  sum  of  squares  (e  e)  relative  to  the  increase  in  regression  sum  of 

aT  A 

squares  (y  y)  which  is  caused  by  adding  (or  subtracting)  new  variables 
0.  in  the  regression.  Ratios  based  on  these  incremental  sum  of  squares  may 
then  be  used  to  determine  the  significance  of  these  added  parameters.  These 
ratios  may  be  formulated  as  test  statistics  for  which,  under  the  assumption  of 
normally  distributed  errors,  e,  standard  significance  testing  may  be  performed. 
That  the  errors  are,  in  fact,  normally  distributed,  is  justified  by  the  central 
limit  theorem  on  the  normality  of  the  distribution  of  a large  number  of  random 
disturbances . 

B.2.1  Significance  of  the  Regression  Equation 

A 

The  solution  0,  Eq.  (B.3),  to  the  minimization  of  error  sum  of  squares 
is  unique  when  X has  full  column  rank.  Any  coefficient  vector  which  differs 

A 

from  the  least  squares  vector  0 must  lead  to  a larger  sum  of  squares.  That 
is,  for  any  other  coefficient  vector,  0, 


Rglpp* 


f 


r 


has  an  error  sum  of  squares 


(y-X0)T(y-X0)  = yT^y  + (0-0)XTX(0-0) 


The  quantity  (0-0)  is  the  "sampling  error"  of  the  regression. 


(B . 13) 


Assume  that  it  is  desired  to  test  the  hypothesis  that  0 is  some  specific 

numerical  value,  say  0q  (e.g.,  0q  = 0 if  there  is  not  dependence  of  the  data 

on  0).  Then,  if  0 is  the  true  value  of  0, 
o 

0-0q  = 0-0  = (XTX)"TXT(X0+e)  - 0 = (XTX)_1XTe  (B.14) 

Using  Eqs.  (B.6)  and  (B.14),  Eq.  (B.13)  then  becomes 


(y-X0o)T(y-X0o)  = eZ//e  + eT(I-^)e 


(B . 15) 


In  order  to  fulfill  the  requirements  for  using  statistical  tests  based  on 

the  sum  of  squares  of  Eq.  (B . 15) , it  is  necessary  to  recall  the  fact  that  an 

indempotent  quadratic  form  in  independent  standardized  normal  variates  is  a 

chi-squared  variate  with  degrees  of  freedom  given  by  the  rank  of  the  quadratic 

form*  (Sec.  15.11,  ref.  [49]).  Furthermore,  if  the  two  normalized  sum  of 

2 2 2 2 

squares  S^/o  and  S^/cr  are  X (v^)  and  X (vj)  distributed,  and  and  S^  are 
independent,  then 


F = 


S./v. 
i i 

S./v. 

J J 


(B . 16) 


is  a Fisher  F ratio  distributed  with  V.  and  V_.  degrees  of  freedom, 


* T 

A sum  of  squares  S.  can  be  written  as  y Ay  where  y is  vector,  A is  a matrix 
of  known  constants.  Then  the  number  of  degrees  of  freedom  of  S.  is  defined 
to  be  the  rank  of  A . 1 
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As  discussed  after  Eq.  (B.6),  the  rank  of  ,Af  is  m-p  and  that  of 
is  p.  Thus,  the  ratio 


eT  (I-c'/y)  s/p 

T 

e e/ m-p 


(B . 17) 


is  distributed  as  F(p,m-p),  if  0 = 0 . But  from  Eqs . (B.6)  and  (B.7), 

T 2 ^ T a rp  A 

e = (m-p)s  and  it  is  easily  shown  that  E (I-v/iOe  = (0-0  )X  X(0-0  ),  if 

o o 

0=0.  Hence, 
o 


(0-0  )X  X(0-0  ) 
o o _ 


F(p,m-p) 


(B . 18a) 


Examination  of  Eq.  (B.12)  shows  that  the  numerator  of  Eq.  (B . 14)  is  the 
difference  between  the  sum  cf  squares  rcgi  easing  on  0q  and  the  sum  of  squares 
of  the  actual  parameters  0.  High  values  of  F(p,m-p)  correspond  to  rejection 
of  the  hypothesis  that  0 = 0 . In  particular,  to  test  whether  the  data  depends 
on  a specific,  parameter,  0^  of  0,  F(p,m-p)  ratio  is  evaluated  for  0q  = 0. 


A recursive  form  of  Eq.  (B.18a)  is 


F(q,m-p)  = 


(0-0  )Xa(0-0  ) 


(B . 18b) 


where  q is  the  number  of  parameters  in  the  regression  at  any  stage  in  adding 
(or  deleting)  parameters  to  improve  fit. 

The  F ratio  tests  are  based  on  quotients  of  regression  "fit"  to  error  "fit" . 
The  quantity. 


„2  _ (X0)A(X0) 
R - T 

y y 

aT  a . T 

= y y/y  y 


(B.19) 


PiOTPMRi 


measures  the  regression  sum  of  squares  to  observation  sum  of  squares.  The 
positive  square  root  of  Eq.  (B . 19) , R,  is  the  multiple  correlation  coefficient. 
The  closer  R is  to  unity,  the  better  the  performance  of  the  subset  of  regres- 
sion variables . Note  that  R is  the  cosine  of  the  angle  between  the  data  vector 
and  the  p dimensional  subspace  spanned  by  the  included  subset  of  regression 
variables . 


From  the  sum  of  squares  decomposition  (B.8),  and  Eq.  (B.19),  it 
may  be  shown  that 


The  closer  R to  unity,  the  stronger  the  dependence  of  the  data  on  the 
regression  parameters  and  the  higher  the  F value.  If  there  is  little  depend- 
ence, R is  "small",  the  hypothesis  Ch  = 0 is  "correct",  and  F is  "small". 

This  test  is  dependent  on  choosing  a critical  value  of  F which  specifies  the 
cutoff  levels.  This  critical  F is  chosen  as  a function  of  m (number  of  obser- 
vations) and  p (number  of  parameters)  from  tables  of  F statistical  tables 
for  a desired  confidence  level.  The  confidence  level  can  be  selected  on  the 
basis  cf  a priori  knowledge  about  level  of  noise  in  data. 


B.2.2  Significance  of  Individual  Conditions  to  the  Regression  Equations 


The  multiple  correlation  coefficient  measures  the  dependence  of  the  data 
on  the  complete  set  of  regression  variables . If  several  sets  of  parameters  are 
to  be  evaluated  with  respect  to  their  respective  explanations  of  the  data,  the 
R coefficient  would  serve  as  one  measure  of  relative  performance.  Unfor- 
tunately, it  follows  from  the  definition  of  R that  increasing  the  number  of 
explanatory  variables  always  increases  R,  unless  the  added  variables  are 
linearly  related  to  other  included  variables.  Beyond  a certain  point,  the 
added  variables  start  fitting  the  random  noise  e.  For  example,  if  the  number 


I 
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of  variables  is  increased  to  m,  a perfect  fit  can  be  obtained  and  R is  one. 
In  order  to  eliminate  parameters  not  significant,  only  the  most  linearly  in- 
dependent variables  are  desired. 


Stated  in  terms  of  the  sum  of  squares  principal,  it  is  desired  to  incor- 
porate those  new  variables  into  the  regression  which  most  reduce  the  error, 
given  that  the  other  variables  are  included  in  the  regression.  If  the  increase 
in  regression  sum  of  squares  due  to  adding  a new  variable,  say  x. , after 

I A 

variables  x.  have  been  included  in  the  regression,  is  denoted  as  Aj  | yj  j , 
then  the  ratio 


t2  = — AlJill 

yxj*x1'X2..-..x._1,xj+1 xp  aTa  + A|1  AJ 


(B . 21) 


measures  the  partial  contribution  of  x.  to  the  regression.  The  square  root 

2 3 

of  R is  the  partial  correlation  coefficient.  Several  forms  of 

yxj  xi’-,,*xp 

this  definition  may  be  written.  One  is  defined  by  letting  y be  the  residuals 
from  the  regression  of  x^  on  the  same  variables.  It  may  then  be  shown 


R 


®«L>‘ 


CKjXEy)2 


(B  .22) 


Note  that.,  in  general,  E R --  f R^. 

j yXj 

Now,  let  the  regression  equation  y = X0  + e be  partitioned  as 


y = X101  ♦ X2e2  + e 


(B . 23) 


where  X^  includes  q variables  and  X2  contains  p-q  variables.  Then 


y'XlGl  = X2°2  + G 


(B.24) 


which  shows  that  an  estimate  of  0^  could  be  obtained  by  regressing  the  residuals 

from  the  regression  of  y on  X,  (i.e.,  which  estimates  0,).  Then  the  vector 

1 ci)  1 

y-X^0j  is  regarded  as  a new  observation,  say  y , which  may  be  regressed 
on  to  estimate  0£  • This  decomposition  can  be  applied  to  each  possible 
subset  of  variables,  X.,  "bringing  in"  new  variables  from  the  right  to  left 
hand  side  of  Eq.  (B.24).  The  requirement  on  "bringing  in"  new  variables 
may  be  satisfied  by  examining  the  significance  of  each  variable. 


The  F test  may  be  used  to  determine  the  significance  of  a single  para~ 

2 2 2 2 

meter  by  noting  the  estimate  of  the  variance  o , s , is  distributed  as  o % 

2 2 2 m p 

Hence,  s /o  ~ Cx  -p^/(  - ) ^rom  Eq-  (B . 7) . Then  for  the  parameter  0., 


0.-0.  (0.-0.  )/ctq. 

ii  l i °i 


S0./O0. 

l l 


(B  .25) 


where  Sq  is  the  standard  error  of  0.,  which,  from  Eq.  (B.8)  is 


sn  - s J s.. 
0.  v n 

l 


(B . 26) 


T — i 

where  s^  is  the  square  root  of  the  i diagonal  term  of  (Xa) 

A 

Since  (0i-0.)/cr^  -*1(0,1),  it  follows  that,  by  definition  of  Student's  t 
distribution  that 


0.-0. 
l l 


(B . 27) 


In  particular,  if  it  is  desired  to  test  the  hypothesis  0^  = 0 (i.e.,  y does  not 
depend  on  0),  the  statistic  t = Qj/sq  is  used.  It  is  shown  in  [49]  that  the 

i 2 

F distribution  with  1 and  (m-p)  degrees  of  freedom  is  equivalent  to  the  t 
distribution  with  rn-p  degrees  of  freedom.  Hence,  the  significance  of  indi- 
vidual regression  coefficients,  0,  is  determined  from  F ratios 
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(B.28) 


F = 02/s2 
1 1 


If  the  ratio  (B.28)  indicates  a variable  is  not  significant,  then  the 
variable  is  deleted.  To  bring  in  another  variable,  the  partial  correlation 
coefficients  of  all  other  parameters  are  examined.  To  form  the  F ratio  for 
these  coefficients,  Eq.  (B.21)  (with  Eq.  (B.26)  may  be  manipulated  to  show 


R-~  = 


<Vse.> 


yXj  (9  /se  r + (m-q) 
J i 


(B . 29) 


where  q is  the  number  of  variables  already  in  the  regression.  The  corres- 
ponding F test  is 


F.  = 
J 


(m-q) 


(B . 30) 


The  variable  (F  ratio  with  1 and  m-q  degrees  of  freedom),  is  calculated  for 
each  of  the  remaining  variables.  The  variable  with  the  highest  value  is  then 
brought  into  the  regression. 


B.2.3  Summary  of  the  Subset  Regression  Method 

The  variables  x^,X2,...,x  are  postulated  as  possible  causative  factors 
in  determining  the  observed  data,  y (Eq.  (B.23)).  The  variables  x,,xOI...x 

12  p 

are  ranked  according  to  the  magnitude  of  their  individual  partial  correlation 
coefficients  with  y (or  F ratio,  Eq.  (B.38)).*  The  variables  with  the  highest 
significance  (F  to  enter)  is  added  to  the  regression.  The  significance  of  the 
added  parameter  is  then  tested  to  determine  if  it  is  above  the  critical  F value 
for  deletion  (Eq.  (B.28).  As  variables  are  entered,  new  F ratios  are  computed 
frem  the  remaining  variables  not  in  the  regression  since  the  degrees  of  free- 


Thc  selection  procedure  described  here  is  also  known  as  univariate  stepwise 
regression,  for  subset  regression  of  order  1. 
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dom  have  been  reduced.  This  process  is  terminated  when  no  new  variables 
satisfy  the  b ratio  required  to  enter  and  when  one  are  to  be  removed. 

At  each  variable  incorporation,  the  multiple  correlation  coefficient  (R, 
Eq.  (B . 19)) , the  equation  F ration  (Eq.  (B.20)),  and  the  standard  error  of 
the  entered  parameter  (Eq.  (B.26))  are  all  included. 

It  is  noted  that,  in  general,  the  particular  variables  finally  selected 
are  not  unique.  Use  of  orthogonal  variables  would  result  in  uniqueness. 

B . 3 COMPUTATIONAL  CONSIDERATIONS 


The  previous  sections  have  demonstrated  that  both  classical  and  sub- 
set regression  parameters  are  obtainable  from  steps  in  the  solution  of  a set 
of  linear  equations  (ref.  Eq . (B.3)).  In  order  to  reinforce  this  connection, 
consider  the  augmented  data-coeffi cient  matrix  (with  rank  p+1) 


A = 


T 

XX 

T 

y x 


T 

y y 


Using  the  inversion  of  Eq.  (B.3),  it  is  found  that  the  final  element 

, A-1  f 
oi  A , a , is 


f 

a = 


1'  T T - 1 T T T T 
y y - yXX(X  X)  X y y y-0JX  XO 


1 

T T 
(1-R  ) y y 


i.e.,  the  multiple  correlation  coefficient  (and  the  sum  of  squares  ||  £||  ) 
are  directly  obtainable.  Similarly,  it  is  found  that  the  diagonal  element 
of  (X  X)  corresponding  to  the  i1*1  coefficient  of  x,  is 


(X^X)"1  = 1 T 

u ~2.-l 


(1  - P.  )x.  x. 

l l i 
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S 


where  is  the  multiple  correlation  coefficient  of  the  regression  of  x,  on 

k 

X1  ,x2 ’ ‘ ‘ ' ,xi-l ,xi+l ' ' ' ’ ’ Xn ’ It  may  be  further  shown  that  the  variance  of 

* r 

0.  is 


A 2 T -1 
var  0.  = cT(X  X)..1  = 
i u 


2 1’ 

(l-Ppxx. 


Computationally,  the  inversion  of  A is  based  on  the  Gauss-Jordan 

th 

pivot  algorithm.  Let  the  k diagonal  element  of  A be  nonzero.  Applying 


a Gauss  Jordan  pivot  on  this  element  is  a new  matrix  whose  ij^1  element  i£ 


a. . , l .e . , 
D 


a. . = 
!J 


Sij  aik  akj/akk 

i ^ k 

j * k 

_aik/akk 

i / k 

j = k 

akj/akk 

i = k 

j i k 

1/akk 

i = k 

j = k 

The  final  result  of  this  inversion  is  the  matrix  B, 


B = 


7'  - 1 

(XXX)  1 


-0 


imu 


The  recursive  algorithm  of  the  Gauss- Jordan  pivot  sweeps  through  the  A 

matrix,  generating  statistical  parameters  which  guide  the  deletion  and  addition 

of  new  variables . Details  of  this  selection  for  the  Gauss- Jorden  pivot  arc 

found  in  ref.  [50]  One  further  parameter  of  this  method  is  the  tolerance.  This 

is  a parameter  of  the  computation  algorithm  itself  and  is  a measure  of  the 

accuracy  of  the  calculations  and  computer.  If  a.,  is  the  value  of  the  itJl 

diagonal  element  of  the  non-negative  matrix  (a.  ) after  pivoting  on  several 

elements,  then  a.. /a.,  is  the  tolerance, 
n n 
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APPENDIX  C 
INPUT  DESIGN 


C.l  INTRODUCTION 

Many  new  developments  are  being  made  in  the  field  of  input  design. 
Recently,  two  efficient  algorithms*  were  developed  for  linear  systems.  These 
techniques  are  capable  of  finding  inputs  which  minimize  a weighted  sum  of 
covariances  of  parameter  estimates  or  any  one  of  many  other  criteria  related 
to  estimation  errors  from  the  resulting  experiment.  They  are  referred  to  as 
Frequency  Domain  and  Time  Domain  methods  and  are  summarized  in  Section 
C.2.  For  details,  see  [ 50  ] . 

The  input  design  techniques  for  lineal-  systems  are  modified  to  work 
for  certain  nonlinear  systems.  However,  a basically  different  approach  is 
required  to  design  optimal  inputs  which  will  identify  the  parameters  of  certain 
nonlinear  systems  accurately.  These  techniques  are  described  in  Sections  C.3 
and  C . 4 . 

C . 2 INPUT  DESIGN  FOR  LINEAR  SYSTEMS 
C.2.1  Frequency  Domain  Method  [51,52] 

Consider  the  state  space  representation  of  a discrete  time  system 

x(k+l)  = 4>x(k)  + Gu(k)  k = 1,2,...,N  (C.l) 

and  the  noisy  measurements 

y(k)  = Hx(k)  + v (1< ) (C.2) 


where 


This  work  was  performed  with  a support  from  the  NASA  Flight  Research 
Center,  Edwards,  California,  under  Contract  NAS  4-2068. 


x(*)  is  a n x 1 state  vector 
u(')  is  { x 1 control  vector 

y(*)  and  v(*)  are  m x 1 measurement  and  noise  vectors, 
respectively . 

tp,  G and  H are  matrices  of  appropriate  dimensions  and  contain  p unknown 
parameters  0.  Fourier  transform  (C.l)  and  (C.2)  to  get 

y(n)  = H(e  - tp)  ^ Gu(n)  + v(n) 

A T(n,0)  u(n)  + v(n)  (C.3) 

As  the  number  of  sample  points  increase,  the  average  information  matrix  per 
sample  approaches 


M = 


2n 


9T* 

90 


S~*  (to) 


9T 

90 


dF  (to) 
uu 


(C.  4) 


where  F is  the  spectral  distribution  function  of  u and  S is  the  spectral 
uu  vv 

density  of  v.  The  following  algorithm  converges  to  the  optimal  input  (with 
energy  constraint  on  input) : 


(a)  Choose  a nondegenerate  input  f (to)  (i.e.,  consisting  of  more  than 
^ frequencies , with  a finite  power  in  each  frequency) . 

(b)  Compute  the  function  vj/(co,f)  and  find  the  value  of  to  where  it  is 
maximum,  call  it  to  , where 

i|/(co,f)  = RejTr  s;]v(cc)  ||  M_1  §§-*]  (C.5) 

to  minimize  | D|  , where  DAM  ' 


and 
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Discretize 


Compute  T(n,G) 
and  (n,6) 


Compute  Information 
Matrix  and 
Dispersion  Matrix 


Find  a.  To 
Maximize  |Kr ^ 

. = oM+(l-a.)M(n.) 


Update  Design 
~ o.f.+(l-a.)f(n.) 


Convergence  Criteria 
Is  CL  Approaching  One 

Is  |M44lHM.|<e? 


Practicalize 


Check  to  Ensure  Max  \J» 
Is  Reasonable. 


SlOP 


Figure  C.l  Optimal  Input  Design  Program  in  Frequency  Domain 

(Maximize  I Ml) 


Input 

F.G.H.At.N.O 
Vv<n).  Sutl("> 

n=-N/2 N/2-1 

1 

' 

>,f)  = aji 


D(f)  ^ls-l((0)  D(f) 


(C.6) 


to  minimize  a linear  functional  R of  D. 


(c)  Evaluate  the  normalized  information  matrix  at  (0 


(d)  Update  the  design 


f.  = (1-a  ) f + a f(to  ) 0 < a < 1 

1 o o o o o 


(C.  7) 


is  chosen  to  minimize  |D(f)|  or  [ .^(DCf^))] 


where 


M(£, ) = (1-  a ) M(f  ) + a M(<o  ),  0 < a < 1 (C.8) 

I o o o o o 


It  can  be  shown  that  such  an  a exists. 

o 


(e)  Repeat  steps  (b)  - (d)  until  desired  accuracy  is  obtained. 


Examples  of  application  of  this  method  to  aircraft  problems  are  given  in 
reference  [ 52  ] . 

C . 2 . 2 Generalized  Time  Domain  Method  [53] 

The  continuous  time  representation  of  a linear  system  is 


x = Fx  + Gu 


x(0)  = 0 


o < t < T 


(C.  9) 


x(t)  is  a nxl  state  vector,  u(t)  is  a qxl  control  vector  and  F and  G are 
matrices  with  appropriate  dimensions  and  depend  on  p unknown  parameters 
0.  There  are  noisy  measurements  of  m linear  combinations  of  state. 


H,  an  raxn  matrix  is  a function  of  0 and  v(t)  is  a white  noise  vector  with  zero 
mean  and  covariance  R.  There  is  an  energy  constraint  on  the  input 


The  information  matrix  for  the  unknown  parameters  0 is  given  as 


We  explain  without  proof  how  different  criterion  can  be  optimized  in  this 
generalized  time-domain  approach.  For  details,  see  [ 53  ] . 


Trace  of  the  Information  Matrix* 


A weighted  trace  of  the  information  matrix  can  be  maximized  by  a simple 
transformation  of  parameters;  sec  ref.  [50  ].  For  multi-input,  multi-output 
systems,  this  problem  was  first  solved  by  Mchra  [55). 


Fo  = 


Thus,  we  want  to  maximize  (C.14)  subject  to  state  equations  (C.15)  and  con- 
straing  (C.ll).  The  necessary  conditions  for  the  optimum  u involve  solving 
the  two-point  boundary  value  problem 


"^GeGel  fx6 


nWHo 


X X 


Xq(0)  = 0 


X(T)  = 0 


“opt  = '^Ge  x <c-2: 

The  smallest  value  of  constant  p,  for  a nontrivial  solution  to  the  two-point 
boundary  value  problem  gives  the  optimal  solution.  Constraint  (C.ll)  is  met 
by  rescaling  x^  and  X This  problem  is  solved  by  eigenvalue-eigenvector 
decomposition  techniques  [54]  . 

Determinant  or  Linear  Functional  of  Dispersion  Matrix 

The  following  algorithm  converges  to  a stationary  point. 


(a)  Choose  any  input  u (t),  which  gives  nonsingular  information 


matrix  M . 

o 
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(b)  Find  an  input  u (t)  to  maximize  <p(u)  where 

m 


cp  = Tr(MQ  M)  to  minimize  the  determinant  of  dispersion 


matrix 


= #'(M  * M M *)  to  minimize  a linear  functional  of  the 
o o 


dispersion  matrix 


(C. 22) 


Both  of  these  criteria  can  be  recast  as  maximizing  a weighted  trace 
of  the  information  matrix. 


(c)  The  information  matrix  for  input  auQ(t)  + Pu^ft)  is 


M,  = a2M  + p2M  + 2apM 
1 o r m r om 


(C.  23) 


where  M is  the  information  matrix  for  input  u (t)  and  M is 
m c m om 

a mixed  information  matrix  defined  in  reference  [ 53  ] . The  energy 
constraint  on  the  input  requires 


a2  + p2  + 2ap 


r uT(o  u 

I O ' T 


(t)  dt  = 1 


(C.  24) 


use  (C.23)  and  (C.24)  to  find  a between  0 and  1 which  optimizes 
the  criterion  function  for  M^.  Such  an  a exists. 


(d)  Find  Mj  and  return  to  (b)  if  the  error  criterion  is  not  met. 


(e)  Check  if  the  solution  is  globally  optimum. 


C.  3 NONLINEAR  SYSTEMS 


Methods  for  designing  inputs  to  identify  parameters  of  nonlinear  dynam- 
ical system  depend  among  other  things  on  kind  of  nonlinearity,  availability  of 
computational  optimization  programs,  and  the  way  the  input  is  applied  during 
the  flight  test.  Many  techniques  lead  to  complicated  algorithms,  which  cannot 
be  solved  for  real  problems  in  a reasonable  computation  tinje.  We  start  with 
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general  nonlinear  problem  and  will  later  discuss  algorithms  which  can  be  used 
with  specific  nonlinearities. 

C.3.1  General  Nonlinear  System 

Consider  a system  with  n states  and  m measurements,  which  follows  the 
equations , 


x = f(x  ,u , 0 ,t) 


0 < t < T 


(C.  25) 


y = h(x,u,t)  + v 


(C. 26) 


v is  defined  as  in  (C.ll).  For  the  sake  of  simplicity,  h is  assumed  independent 
of  0.  The  information  matrix  for  parameters  0 is 


■n 


dh(x,u,t) 


-1  /dh(x,u  ,t)\ 

\ d0  I 


(C .27) 


=/T(H)T(i)Vli 


(C .28) 


The  equation  governing  the  sensitivity  of  state  to  parameters  is 


d 9x  _ 9f  . 9f  Ox 
dt  90  90  9x  90 


(C.29) 


We,  thus,  wish  to  choose  a bounded  energy  input  to  optimize  a certain  function 
of  M of  (C.28)  subject  to  the  state  equations  (C.25)  and  sensitivity  equations 
(C.29).  To  maximize  the  trace  of  the  information  matrix,  it  is  necessary  to 
solve  the  following  two  point  boundary  value  problem. 
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xQ  (0)  = 0 


where  x0  is  defined  by  (C.16)  and 


3h  T -1  8h 
3x  K 3x 


f(x,u,0 

,t) 

ar  df 

9x 

30  x 3x 

-|_ 

1 

391 

i 

i 

l 

3 f 3f 

3x 

30  3x 

30 

P 

P 

0 

ah  T -i  ah 

3x  3x 


(C.  31) 


(C . 32) 


The.  optimal  input  is. 


(C.  33) 


(C.30)  is  an  extremely  difficult  problem  to  solve  computationally.  An  approxi- 
mation is  to  find  states  and  controls  to  optimize  the  desired  function  of  the 
information  matrix,  while  satisfying  the  sensitivity  equations  only  approximately, 
The  e-technique  is  a possible  method.  In  this  method,  the  following  function  is 
optimized . 


(C . 34) 


min  J 
x,u 


j j-xeRe  xg  + eilxe  fe^xe'u,0,t^lw|  dt 


The  constraint  accruing  from  the  sensitivity  equations  can  be  tightened  by  in- 
creasing e.  Gradient  methods  can  be  used  to  optimize  J. 

C.3.2  Application  to  Aircraft  Problems 


In  flight  tests  of  most  aircraft,  it  is  customary  to  measure  all  states  and 

accelerations , If  x is  noisy  measurements  of  state,  with  noise  v,  and  a is 
m 1 m 

the  noisy  measurements  of  accelerations  with  noise  v^, 


or 


am‘V2  = f(^m"vl>u'0’t)  CC.35) 

am-v2  = (x,u  ,9  ,t)A0  + (xm,u,0,t)v1  + f(xm,u,0,t)  (C . 36) 

•if  »v/ 

am"f(xm«u>0't)  = ae  (x,u,0,t)A0  + v2+  ^ (xm,u,0,t)v1  (C.37) 


Assuming 

fe(xm,u,e,t)  = H (x,u,0,t)  (C.  38) 

(C.37)  becomes 

y = ll  (x,u,G,t)A0  + v2  + U,G,u,t)v1 

A h(x  ,u  ,0  ,t)A0  + g(x,u,0,t)v  (C.39) 


The  information  matrix  for  parameters  0 (or  AO)  is 


M - (gRwgr)  *h  dt 


(C . 40) 


The  following  procedure  can  be  used  to  determine  an  input  to  approximately 
minimize  the  determinant  or  a linear  functional  of  the  dispersion  matrix. 

(a)  Find  functions  h and  g of  (C.39). 

(b)  Choose  an  energy  bounded  input  history  UQ(t)  and  state  history 

x (t)  and  find  the  corresponding  information  matrix  M . 
o o 

(c)  Find  input  u(t)  and  state  x(t) , under  same  energy  constraints 
as  above  to  maximize  <p(u,x)  defined  in  (C.22). 


(d)  Find  a,  optimal  input  and  states  as  in  step  (c)  of  the  algorithm 
for  linear  system  (see  page  199)  . 

(e)  Find  M^1  and  return  to  (b)  if  the  error  criterion  is  not  met. 

(f)  x°(t)  and  u°(t)  obtained  this  way  will,  in  general,  be  incompatible. 
Find  compatible  input  and  state  by  minimizing 


x°(t)llw  + llu^)  - u°(t)Hw  | dt 

1 2) 


(C.41) 


Subject  to  state  equation  constraints  (C.25).  This  is  a fairly 
simple  optimization  problem. 

C.4  NONLINEAR  SYSTEMS  WITH  SPECIFIC  NONLINEARITIES 


Some  input  design  techniques  can  be  used  only  with  a particular  class 
of  nonlinearities.  Two  such  techniques  are  presented  here. 


C , 4 . 1 Describing  Functions  Method 


This  technique  is  a unique  technique  in  that  the  nonlinear  system  is 
approximated  by  a linear  system  for  design  of  input.  This  method  works 
successfully  only  with  very  simple  nonlinearities  in  low  order  systems.  The 
inlinear  system  is  converted  into  its  "frequency-domain"  representation  by 
one  of  the  describing  functions  methods.  The  method  described  in  Section 
C.2.1  is  then  used  to  determine  the  spectrum  of  optimal  input  in  steady  state. 

C . 4 . 2 Series  Expansion  of  Nonlinear  F u notions 

Every  nonlinear  function  can  be  expanded  in  an  infinite  series  of  a set 
of  complete  functions.  The  truncated  series  is,  usually,  a good  approximation 
for  a finite  range  of  the  independent  variables.  In  parameter  identification 
technology,  it  is  customary  to  express  a nonlinear  function  as  a polynomial  in 
independent  variables  when  the  nature  of  nonlinearity  is  unknown  or  difficult 
to  describe  exactly.  In  the  high  an gle-of- attack  flight  regime,  the  aerodynamic 
derivatives  are  functions  of  angle-of- attack , Mach  number,  altitude,  etc. 

Each  stability  and  control  derivative  is  expanded  into  a suitable  order  poly- 
nomial in  the  independent  variables  (say  angle-of-attack , a,  and  Mach  number 
M) , with  unknown  coefficients,  viz. 

C (a,M)  = P (M)  + P,  (M)a  + P_(M)a2,  . . . ,P  (M)ak  (C.42) 

m o x & k 

q 

and  P^(M)  are  suitable  polynomials  in  mach  number. 

P.  (M)  = C1  + clM  + C^M2  + ...  + CoM£  (C.  43) 

i o 1 2 £ 

The  problem  of  finding  the  nonlinear  function  of  C is  reduced  to  deter- 

m 

q 

mination  of  the  values  of  the  polynomial  coefficients.  The  order  of  the  poly- 
nomial depends  on  the  kind  of  nonlinearity  and  the  range  of  the  independent 
variables  for  which  the  function  is  to  be  defined. 

I 
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Each  different  set  of  complete  functions  lead  to  a slightly  different 
method  for  determining  inputs  which  will  accurately  identify  the  coefficients 
of  the  terms  in  the  expansion  accurately.  We  discuss  here  various  aspects 
of  input  design  for  systems  which  are  linear  except  that  some  of  the  coef- 
ficients are  expressed  as  polynomial  expansions  in  some  state  variables. 


The  ultimate  objective  of  any  input  is  to  determine  the  parameters  of 
the  expansion  such  that  a good  definition  of  the  function  is  obtained  over  the 

range  of  interest.  Let  us  see  how  it  can  be  achieved.  If  & be  the  error  in 

i . 1 

estimating  , the  error  in  the  value  of  the  aerodynamic  coefficient  would  be 


(for  example  in  C ) 

mq 


and 


Cm  (ct,M)  = Pq(M)  + P (M)u  + ...  + F,  (Ivl ) a 
mq 


P.(M)  = C1  + C\M  + ...  + Cl„M* 
1 o 1 x 


(C . 44) 


(C . 45) 


K 


Thus,  to  get  minimum  mean  square  error  with  weight  d(a,M)  in  a certain 
range  £2  of  a and  M,  it  is  necessary  to  minimize 


J =J  d(a,M)  jcm  (a 
(a,M)  efi 


,M)j‘ 


dordl\/l 


i 


= Tr(WM-1) 


(C .46) 


where  M is  the  information  matrix  for  the  coefficients  C?  (order  C°,C? 
o 1 I ] k J ° 1 

/-.i  /~-h\  _ i 


' ’ ’ o'  o'  i,<22  ’ ’ ' ’ *Cj) * and 


W 

ki+j.ki'+j 


«+1  '2Mj+j  "2d(«,M)da  dM 
(a,M)t;fl 


(C. 47) 
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There  are  two  basic  methods  for  identifying  the  coefficients  of  polynomial 
expansions.  In  the  first  method,  these  polynomial  coefficients  are  extracted  from 
one  input.  This  requires  inputs  with  suitable  amplitudes  so  that  the  independent 
variables  of  the  polynomial  expansions  vary  throughout  the  region  of  interest. 
The  other  important  consideration  in  this  method  concerns  the  order  of  the  poly- 
nomial which  can  be  identified  with  enough  confidence.  This  method  is  good 
where  the  variation  in  independent  variable  is  small  and  can  be  achieved  within 
a reasonable  time. 


The  second  method  identified  the  value  (and  maybe  slope)  of  the  function 
at  different  values  of  the  independent  variable.  This  is  done  by  exciting  the 
aircraft  modes  at  the  desired  trim  value  of  the  independent  variable.  A poly- 
nomial can  then  be  passed  through  these  points  by  a least  squares  type  tech- 
nique. Here  it  is  necessary  to  select  the  values  of  the  independent  variables 
at  which  the  experiment  should  be  carried  out  and  also  the  duration  of  the 
experiment  at  each  point. 


C.4.2.1  Input  Amplitudes  to  Identify  Polynomial  Expansions 


In  linear  systems,  parameter  estimate  error  covariances  scales  with 
input  energy.  If  the  energy  is  doubled,  the  variance  in  all  parameter  esti- 
mates decreases  by  a factor  of  two.  In  nonlinear  systems,  this  is  not  true. 

In  particular,  the  accuracy  of  the  coefficients  of  high  order  terms  in  the 
expansion  increases  in  much  higher  proportion  with  increasing  state  and  input 
amplitudes.  For  very  small  amplitudes,  the  system  "looks"  like  a linear 
system  obtained  by  linearizing  the  nonlinear  system  around  the  average  state. 


To  simplify  the  discussion,  we  assume  that  all  states  and  acceleration 
are  measured  and  that  the  system  involves  only  one  nonlinear  derivative  C 
which  is  expressed  as  a polynomial  in  some  state  variable  a.  If  the  measure 
ments  are  substituted  in  state  equation  (C.25) 


where  is  the  measurement  noise  in  x,  etc.  The  above  equation  can  be 
solved  for  C(a)  to  give 


let 


y (t)  = C(a(t))  + n(t) 


C(a)  = C..  + C_  a + ...  + C a™'1 
i c.  m 

and 


( 


( 


An  estimate  of  C from  this  experiment  is. 


where  r is  the  power  spectral  density  of  n.  The  covariance  of  C is  lower 
bounded  by 


Cov(C)  > 


M 


■/.i 


■/ 


1 

a - • 

m-1 

• • • a 

a 

2 

m 

a 

a 

m-1 


da 


b d(a) 


1 a 


a a 


m-1 


2m~2 


m-1 


m 


a2m-2 


da 


(C.  54) 


(C.  55) 


where  a varies  between  a and  b and  d(a)  is 


d(a*)  = I 

whenever  a 
a=a* 


(C . 56) 


d(a)  is  a measure  of  the  fraction  of  time  spent  in  the  neighborhood  of  a 
during  the  experiment.  To  study  the  nature  of  estimation  error,  it  is  best 
to  expand  C(a)  in  polynomials  which  are  orthogonal  over  the  interval  a < a < b 
with  kernel  d(a),  i.e.  , 


C(a)  = C'P.  (a)  + ClP-(a)  + ...  + C'P  (a) 

i i c.  c.  mm 


CC.57) 


where  P.(a)  is  a (i-l)st  order  polynomial  in  a.  Then  the  information  matrix 
for  parameter  vector  C'  is 


M'  = 


d(o) 


Pl(c0 


Pj(a)P2(a) 


P (a)P  (a) 
l m 


Pj(a)P2(a) 

P2(«) 


P1  (a)pm  (a) 


Pm(a) 


da 


(C.  58) 
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P?(a)da 


Since  from  condition  of  orthogonality 


P.  (a)  P. (a)da  = 0 


P (a)da 


Thus,  the  variance  of  C 


d(a)pf(a)da 


cov(C.')  > 


The  ratio  of  parameter  value  to  its  standard  deviation  is 


b d(a)P.  (a)  da 


A good  d(a)  can  be  established  based  on  the  system  and  the  class  of  inputs. 
Then  the  t-ratio  for  each  parameter  is  determined  for  the  experiment  duration 
The  t ratios  usually  decrease  with  increasing  i and  decreasing  (b-a)  . A 
suitable  (b-a)  is  chosen  to  obtain  acceptable  t-values  for  the  coefficients  of 
the  polynomial.  For  the  case  where  equal  time  is  spent  at  all  a (a  < a < b). 


and  P^(a)  are  (for  a = -1,  b = 1), 


Po(a)  = 1 


P^a)  = x 


P2(a)  = x6  - 1/3 


P3(a)  = x3  + 5/3x2  - 5/9 


The  expression  on  the  RHS  of  (C.62)  can  be  computed  for  all  i with  r as  a 
parameter . 


C.4.2.2  Selection  of  Experiment  Points  to  Identify  Polynomial  Expansions 


Consider  the  system 


x ~ F(x)  x + G(x)  u 


1 


The  matrices  F and  G contain  aerodynamic  derivatives  which  are  expressed  as 
polynomial  expansions  in  some  state  variables.  Suppose  an  experiment  is 
carried  out  starting  from  a trim  value  xq.  It  is  assumed  that  the  terms  in 
F(x)  and  G(x)  are  slowly  varying  functions  of  x,  so  that  during  this  experi- 
ment these  matrices  are  almost  constant.  We  can  design  optimal  input  for 
this  linear  systems  so  as  to  obtain  "best  estimates"  of  F(x)  and  G(x). 

The  accuracy  of  these  estimates  depends  on  xq  and  the  time  t spent  during 
the  experimentation.  In  particular  for  an  aerodynamic  coefficient  C which 
is  a function  of  some  state,  say  Q, 


C(a  ) = C,  + C.,cr  + . 

o e 1 2 o 


, „ m-1  . . 

. + C a + n(a  ,t  ) 
mo  o o 


where  the  noise  nCa^.t^)  has  mean  zero  and  covariance  which  depends  on 

a and  decreases  linearly  with  t , i.e., 
o 7 o 
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*<»  <w> 


t 

o 


( 


It  is  clear  that  if  the  experiment  is  carried  out  in  the  neighborhood  of  x. 
(i=l , 2 , 3 , . . . E ) for  time  t. , the  coefficients  Ch's  can  be  estimated  using  the 
formulae 


A 

C A 


= M_1Y 


( 


where 


The  Cramer-Rao  lower  bound  on  parameter  estimate  covariance  is 


cov(C)  > M 


The  covariance  of  the  estimates  approached  M ^ asymptotically.  To  obtain 


good  estimates  of  the  nonlinear  function  over  the  range  of  interest,  it  is 
necessary  to  minimize  a weighted  trace  of  M 1 with  the  constraint  of  fixed 


total  experiment  time , 


£ t.  = T 
i=l  1 


The  information  matrix  has  the  following  properties: 


(a)  It  is  a positive,  semi-definite  matrix  as  is  clear  from  (C.70). 


(b)  The  information  matrix  is  sum  of  fi  matrices  each  of  rank  one. 

Therefore , the  rank  of  M cannot  exceed  £ . In  other  words , the 
information  matrix  is  singular  if  the  number  of  experiment  points 


is  fewer  than  m. 


We  now  state  and  prove  a theorem  which  leads  to  an  algorithm  for  choosing 
fi,  oc  and  h (i=l ,2 ,3 , . . . ,E) . The  set  of  £ , Ct,  and  t.  is  called  a design 


Theorem  C.l.  The  following  are  equivalent: 


(a)  A design  a*  minimizes  Tr(WM*  ) 


(b)  A design  «*  minimizes  maxfl  { cp  * (a) } 


(c.)  A design  a*  is  such  that 


maxQ  cp*(a)  = Tr(WM*  ) 


where 
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<p(a)  = T^WM_1 


1 

a • • • 

■ • am 

2 

m 

a 

a 

a 

m-1 


a 


2m- 2 


m 


-1 


(C.  75) 


Proof:  Consider  a design  a*  which  minimizes  Tr(WM*_1).  The  information 

matrix  for  a design  which  a linear  combination  of  L*  and  another  design  in 
which  time  T is  spent  at  a is 


M = 2) 
i=l 


a. 

i 


m-1 


a. 

i 


t.(l-p) 

r(u.) 


(1,  a., 


m-1- 
. a.  ) 


r ~i 
l 

a 


m-1 


TCP) 

r(a) 


(1,  a,  • 


m-lx 

, a ) 


(C.  76) 


Since  L*  minimizes  Tr(WM  ^) 


d Tr(WM  *) 


ap 


p = 0 


- Tr(WM  1 M'1) 

dft 


> 0 


P = 0 


(C.77) 


Also 
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-V-.  (1,  a., 

rCap  i’ 


m-lx 
• ■ Bi  ) 


om'1 

- 1 — 


+ T 
rlaT 


2m-  2 


From  (C. 77)  and  (C.7S) 


+ TrCWM*"1) Xr-  Tr  WM*  1 la  am  1 M 1 > 0 

r fa;  — 


2m~2 


cp*(a)  < TrCWM11”1) 


Since  we  can  choose  any  a from  the  range  of  interest,  (C.79)  holds  for  all  a. 
Also  for  any  design  L 


£ cp  (a)  -ip 


, E t.  1 a 
= Tr  WM  I -iv 

a=0l^  a a2 


am  1 a2m"2 


= Tr(WM_1MM  !)  = TrfWM-1) 
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Since 


£ t. 

Z and  t.  > 0 (C.81) 

i=l  X 1 


max^  cp  (a)  > Tr.CWM 


(C.  82) 


(C.79)  and  (C.80)  are  in  contradiction  unless 

maxa  cp(a)  = Tr(WM*_1)  (C.83) 

which  proves  (c)  and  in  the  light  of  (C.82),  proves  (b) . The  assertion  that 

(b)  or  (c)  imply  (a)  is  quite  straightforward. 

The  following  procedure  can  be  used  to  find  the  optimal  design  L*. 

(a)  Evaluate  r(a)  based  on  design  of  linear  experiments  in  the  neighbor 
hood  of  a.  This  function  may  also  include  "cost"  of  carrying  out 
the  experiment  about  state  a. 

(b)  Choose  a nondegencrate  design  L.  It  should  have  at  least  m ex- 
periment points  where  finite  time  is  spent.  Find  M ^ . 

(c)  Find  or  where  cp(a)  is  maximum.  This  cp  (a)  is  more  than  Tr(YVM 

(d)  Find  a p between  0 and  1 such  that  by  spending  time  (T1~P)  in 
points  of  design  L and  Tp  at  a,  the  Tr(WM  ^ ) is  minimized. 

(e)  The  new  design  has  £ + 1 points  at  or.  (i=  1,2 fi),  or  with  time 

(1~P)L  (i=l ,2,  ...,£) , Tp  at  these  points  respectively. 

(f)  Check  if  max^  <p  (o)  is  close  to  Tr(WM  ^ ) . If  not,  return  to  (c). 


The  fact  that  max^  cp(a)  is  greater  than  Tr(WM  is  clear  from  (C.82). 

Also,  from  (C.77)  and  (C.78),  for  small  p 

Tr(WM_1)p  = Tr (WM_1)p_0  - p(cp(a)  - Tr (WM-1) p=())  (C.84)  ! 

Since  Tr(WM  decreases  for  small  p,  it  is  clear  that  p of  step  (d)  exists. 

Also,  as  long  as  the  design  is  not  optimal,  it  is  always  possible  to  decrease 
Tr(WM  ^).  Since  it  is  bounded  from  below,  it  must  converge  to  the  optimum. 

It  can  also  be  shown  that  the  number  of  points  in  the  design  cannot 

, m(m  + l)  T . . 

exceed  ~ . In  many  cases,  it  is  m. 


APPENDIX  D 

SIMULATION  EQUATIONS 
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This  appendix  is  drawn  from  reference  [1]  wherein  the  complete  simu- 
lation  used  for  generating  the  data  is  described.  Because  of  the  significance 
oi  the  detail  of  these  equations,  and  the  need  to  refer  to  them  in  this  report, 
this  appendix  is  given . 

A flow  chart  of  the  simulation  is  given  in  Figure  D.l.  The  following 
features  arc  notable: 

1.  The  aerodynamic  data  are  "table  look-up"  and  valid  for 

-10  < a < 110  and  -40  < p < +40°.  These  data  are  reported  in 

reference  [ ] . 

2.  Three  types  of  controls  arc  available— pilot  stick  (autopilot), 
pilot  stick  (read  in)  and  stability  augmentation  system. 

3.  The  time  histories  of  responses  and  aerodynamic  coefficients  may 
be  outputted.  The  time  histories  are  either  outputted  as  calcu- 
lated (ideal)  or  they  are  passed  through  nonlinear,  randomly 
disturbed,  instruments. 

Further  program  characteristics  are  listed  in  Table  2.1. 

D.l  EQUATIONS  OF  MOTION  IN  BODY  AXES 
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Figure  D.l  Flow  Chart  of  Six  Deg ree-of- Freedom  Simulation 
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D.1.2  Rotational  Equations 
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D.1.3  Stability  Augmentation  System 


8s1  = .15q  - 8s1 


pitch 


6r'  = 3r  - ^ ' 
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if  |8r'|>5  , set  6r'  = (sgn  6r')5 


Sr'  = 8r  ' + ,0168a  g 
1 yb 

6a1  = ,265p  if  |6a’|  >7.5°. 

if | 6s ' | >0?  5,  set  6s'  = (sgn  6s ’ ) 0?  5 
D.1.4  Autopilot 

v 

8a  = 7r  + cpn 
r 

6a  = -3(PB 
ca 

8 = + v 

s ' 

D.1.5  Total  Controls 


where 
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6!-’  6I»  are  input  controls;  8a,  8a,  8 

X a b 7T  * Si  * 

8' , 8'  8'  are  SAS  controls;  8 ,8  ,8 


set  8a'  = (sgn  6a')7.5 


< 30° 

< 30° 

< 8s  < 9° 

a are  autopilot  controls; 
are  trim  controls. 


D . 2 INSTRUMENTATION  EQUATIONS 


D.2.1  Attitude  Gyros 


scale  factor  on  pitch 


D.2.2.1  Resolution 


where  vl/  , 0 . v]/ 
P P 

misalignment . 


9 are  random  variables  of  gyro  Euler  angle 


;urement: 


D.2.3  Angular  Accelerometers 


D.2.3.1  Resolution 


D.3.3.2  Measurements 


D.2.4  Velocit 


where 


M = Mach  no 


D.2.5  Angle-of- Attack 


D.2.5.1  Indicated 


where 


are  random  variable 


where 


is  a constant,  estimated  distance  from  the 


1 


D.2.6  Sideslip  An  vie 


D.2.6.1  Indicated 


, V sin  p + r£  - p£ 

P.  = (1  + kft)  tan  1 r~ - - + b + w 

1 R Vg  cos  a cos  p+q£^-r£  P ^ p 


z y 


D.2.6. 2 Corrected 


r £ 

p = p.|cos  a |-  — :t7— • cos  p. 
c i'  c'  V . 1 i 


D.2.7  Linear  Accelerometer 


D.2.7.1  CG  Acceleration 


u - vr  + wc 


v - wp  + ur 


w - uq  + vr 


D.2.7. 2 Accelerometer  Triad 
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where  x , y , . . . , z are  random  variables 
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D.2.7.3  Axial  Accelerometer 
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where 


^A’  are  ranclom  variables;  Euler  angles  of  accelerometer 
misalignment 


y * Zr-T  are  constants,  estimated  position  of  the  accelerometer 
CgA  CgA 


D.2.7.4  Lateral  Accelerometer 
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where  cp^,  are  random  variables,  Euler  angles  of  accelerometer 

misalignment 


X z 

cgj  ’ ‘cgj.  are  constants,  estimated  position  of  the  accelerometer 


— 4 
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D.2.7.5  Vertical  Accelerometer 


where 


are  random  variables,  Euler  angles  of  accelerometer 
misalignment 

are  constants,  estimates  of  the  accelerometer  position 


D.2.7.6  Pilot  Station  Vertical  Accelerometer 


where 


Xps  ’ ^ps  * Zps  are  ranc^oin  variables,  the  location  of  the  pilot 
station  relative  to  the  c.n  . 
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APPENDIX  E 

IDENTIFICATION  MODEL  EQUATIONS 


E.l  INTRODUCTION 


The  equations  of  motion  for  the  aircraft  simulation  are  given  in  the  pre- 
ceding Appendix  D.  These  equations  are  to  provide  the  framework  for  data 
generation.  For  parameter  identification,  a less  general,  more  approximated 
model  representation  is  desired,  which  still  retains  the  essential  features  of 
the  known  aircraft  dynamics . Similar  approximations  must  be  applied  to  the 
instrumentation  models.  This  appendix  outlines  the  models  used  in  the  iden- 
tification program. 


E . 2 DYNAMIC  AND  MEASUREMENT  EQUATIONS 


The  state  equations  are: 

qBS 


— A A A A 

u = vr  - wq  + 
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C + C 8 + — — C q + — 

x Xo  s x M m 

o 2V~  q 

s B H 
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..  A A A A 
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q*s 


w = uq  - vp 


A A -B" 
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C+C.  8 + C q 
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+ + g cos  0 sin  (p 
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(E . 3) 


and,  for 


C£  - C£  + C 8r 
o 
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+ 6 +4-  (c  r + C p\ 

6a  2VB  \ r P / 


(E.4) 


(E.5) 
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: ] 


sin  <p  + r cos 


The  measurement  equations  are 


Comparison  of  Eqs . (E.13)  - (E.12)  with  the  detailed  nonlinear  instrumentation 
models  of  Appendix  D indicates  the  degree  of  approximation  of  the  measurements 


E.  3 SENSITIVITY  EQUATIONS 


(E.ll),  may  be  functionally  written 


The  dynamic  equations 


where  C is  a coefficient  matrix  of  constants  (which  may  or  may  not  be  known), 
f are  inertial,  engine  gyroscopic  and  gravitational  terms,  f are  engine  thurst 
terms,  and  f&  are  aerodynamic  terms.  Denoted  9.  as  a particular  parameter, 
the  sensitivity  of  Eq.  (E.21)  to  0.  is 


PC 

00. 


(f 


+ f + f ) 
a e 


where  Of  /00.  is,  in  general,  zero.  Note  that  the  assumption  of  Of  /O0.  = 0 would  nol 

Cl  G 1 

be  used  for  cases  where  the  engine  thrust  is  highly  dependent  on  state. 

If  the  identification  requirements  include  separation  of  C , C C 

1 1 x«’  y*  z* 

a * p a 

OC/Oth  will  be  non-zero,  otherw.sc  0C/3Q.  = 0. 

The  expression  for  Of^/90^  is 
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v j.  i fi 

vae.  2 0G. 
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where  C are  the  static  coefficients  and  C,^  the  dynamic  rate  coefficients  witJi 
general  expansions 


C = C (a, pi  + £ C + I C (.  p(j)  + HC 

° i cr  ; j p 


i j a(i)p^ 


a(i)p(i) 


and  the  coefficients  0 are  the  coefficients  of  a~'\  p^\  or  cn'^p^. 

Sensitivity  equations  for  the  measurements  are  obtained  from  Eqs . (E.13) 
(E. 20) . 
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E.4  ALTERNATE  IDENTIFICATION  EQUATIONS  FOR  THE  (a,  p,  V)  FRAME 


TABLE  E . 1 

Equations  of  Translational  Motion  in  a,  p,  V System 


TABLE  E . 2 

Accelerometer  Equations  in  a,  p,  V System 


A 

x 

c-g 


= V cos 


P cos  a - V [cp  cos  a + r)  sin  p + (a  - q)  sin  a cos  pj 


A 


y. 


= V sin 


•g' 


P + V jcos  p(P  + r cos  a - p sin  a)  + g cos  0 sin  cpj 


- g sin  0 


A 


= V cos  p sin  tt+  V (p-p  sin  a)sin  p + (ci-q)cos  p cos  a - g cos  O cos  <1> 
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